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not  "sufficiently  long"  to  ensure  that  the  estimation  error  is  small,  estimator 
performance  is  dominated  by  large  estimation  errors  due  to  anomalous  estimates, 
and  actual  performance  can  be  much  worse  than  predicted  by  the  CRLB.  This  report 
presents  a  detailed  analysis  of  time  delay  estimator  performance. 

By  way  of  background,  the  generalized  cross  correlation  (GCC)  approach  for 
TDE  is  discussed  and  it  is  shown  that  many  of  the  commonly  used  TDE  methods, 
including  the  ML  estimator,  are  related  through  the  GCC  method.  In  recent  work, 
a  correlator  performance  estimate  (CPE)  and  a  modified  Ziv-Zakai  lower  bound 
(ZZLB)  have  been  developed,  which  predict  performance  more  accurately  than  the 
CRLB  for  large  estimation  error  conditions.  Derivations  of  the  CRLB,  the  CPE 
and  the  ZZLB  are  given,  and  the  behavior  of  these  performance  estimates  is 
investigated  for  the  case  of  flat,  low-pass  signal  and  noise  power  spectra. 

The  CPE  and  the  ZZLB  are  seen  to  be  characterized  by  a  threshold  signal -to-noi se 
ratio  (SNR).  Above  this  threshold  SNR,  the  CPE  and  the  ZZLB  coincide  with  the 
CRLB,  while  below  the  threshold,  the  CPE  and  the  ZZLB  deviate  from  the  CRLB  and 
predict  much  poorer  performance  than  the  CRLB.  Further,  the  threshold  SNR  is 
shown  to  be  approximately  inversely  proportional  to  the  square  root  of  the 
coherent  processing  time. 

This  behavior  has  significant  implications  for  TDE  signal  processing  tech¬ 
niques.  In  particular,  significant  performance  gains  can  be  realized  by  imple¬ 
menting  a  coherent  processor  as  opposed  to  an  incoherent  processing  algorithm. 
Simulation  results  are  presented  to  corroborate  these  performance  predictions. 

For  the  case  of  a  time-varying  time  delay,  these  observations  point  to  the 
necessity  of  pre-processing  the  received  signals  to  compensate  for  the  relative 
time  compression.  This  compensation  is  required  in  order  to  implement  a  coherent 
processor  and  results  in  a  more  complex  structure  for  the  time  delay  estimator. 

A  relatively  simple  compensation  technique  is  described,  and  preliminary  simula¬ 
tion  results  are  presented  to  demonstrate  the  effectiveness  of  this  technique. 
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CHAPTER  1 


INTRODUCTION 

Time  delay  estimation  (TOE)  is  an  area  of  active  research  with 
applications  in  a  wide  variety  of  fields  [1].  An  important  concern  in  any 
estimation  problem  is  that  of  predicting  the  performance  of  the  estimation 
methods  employed.  It  is  often  possible  to  set  a  measure  of  performance  by 
establishing  bounds  on  performance.  While  performance  bounds  set  an 
absolute  standard  or  limit  for  performance,  the  usefulness  of  a  given  bound 
for  predicting  performance  depends  upon  the  tightness  of  the  bound,  i.e., 
how  close  the  bound  is  to  the  greatest  lower  (or  least  upper)  bound.  This 
research  conducts  a  comparative  study  of  performance  bounds  for  the  time 
delay  estimation  problem.  For  the  TDE  problem,  lower  bounds  on  the 
variance  of  the  time  delay  estimate  can  be  computed.  An  important 
objective  of  this  study  is  to  gain  a  better  understanding  of  the  behavior 
of  these  bounds  as  functions  of  the  signal-to-noise  ratio,  bandwidth,  and 
observation  time.  In  addition,  the  TDE  performance  bounds  are  compared  to 
simulation  results  to  determine  their  relative  merit  for  predicting  time 
delay  estimator  performance.  The  results  of  this  study  are  found  to  have 
significant  implications  for  TDE  signal  processing  techniques. 

1.1  TDE  Background 

In  the  basic  time  delay  estimation  problem,  a  signal  and  a  delayed 
version  of  the  signal,  both  of  which  may  be  corrupted  by  noise,  are 
received  at  a  sensor  or  sensors  and  the  object  is  to  estimate  the  delay 
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value.  This  problem  has  found  applications  in  seismic  prospecting  [2], 
bio-medical  engineering  [3],  nuclear  reactor  engineering  [4],  and  many 
other  areas.  Depending  on  the  application,  the  signal  may  be  electrical, 
acoustical,  or  even  neutron  flux  fluctuations  [4],  and  the  receiver  may  be 
a  single  sensor,  multiple  sensors,  or  sensor  arrays.  The  application  which 
motivates  this  research  is  that  of  passive  sonar  source  localization.  The 
problem  will  be  modelled  as  depicted  in  Figure  1.1.;  an  acoustic  signal 
emanates  from  a  point  source,  propagates  through  the  ocean  medium,  and  is 
received  at  two  spatially  separated  sensors.  Assuming  a  constant 
propagation  speed  of  sound  in  the  ocean,  the  time  delay  in  the  time  of 
arrival  of  the  signal  at  the  two  sensors  is  given  by  the  difference  in  the 
path  lengths  between  the  source  and  the  two  sensors,  divided  by  the  speed 
of  sound  in  the  ocean. 


For  a  given  value  of  the  time  delay,  it  can  be  shown  that  the  source 
must  be  located  along  a  hyperbola  as  shown  in  Figure  1.2.  Taking  the 
baseline  between  the  two  sensors  as  the  x-axis,  with  the  origin  at  the 
mid-point  between  the  sensors,  the  equation  of  the  hyperbola  is  given  by 
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-4-7  =  1/4  ,  (1-1) 

aR^-L"^ 


where  L  is  the  distance  between  sensors  and  aR  is  the  difference  in  the 
path  lengths  between  the  source  and  the  two  sensors.  Note  that  aR  =  D  c, 
where  D  is  the  time  delay  of  the  signal  at  sensor  1  relative  to  sensor  2 
and  c  is  the  speed  of  sound  in  the  ocean. 
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Figure  l.l  Physical  Model  of  the  TDE  Problem 
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Figure  1.2  Bearing  Angie  Interpretation 
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The  bearing  of  the  source,  relative  to  the  midpoint  of  the  baseline, 
is  given  to  a  good  approximation  by  the  bearing  angle,  ©,  defined  by  the 
intersection  of  the  hyperbolic  asymptote  with  the  baseline.  It  can  be 
shown  that 


cos  e 


L 


1 


(1-2) 


where  R  is  the  distance  to  the  source  along  the  hyperbolic  asymptote. 
For  a  distant  point  source  (R  >>  L),  e  is  a  very  good  approximation  to 
the  true  bearing  and  (1-2)  becomes 


e 


cos 


-1  ^ 
L 


cos 


-1  D  c 
L 


(1-3) 


Thus,  if  L  and  c  are  known  and  the  time  delay  between  the  two  sensors  can 
be  estimated,  an  estimate  of  the  source  bearing,  relative  to  the  position 
of  the  sensors,  can  be  obtained  from  equation  (1-3).  For  D  >  o,  the  source 
is  closer  to  sensor  2  than  to  sensor  1  as  in  Figure  1.2.  For  D  <  o,  the 
situation  is  reversed  and  the  source  is  located  along  the  corresponding 
hyperbola  reflected  through  the  y-axis.  The  maximum  possible  delay  occurs 
when  the  source  is  in  the  endfire  position  (col inear  with  the  sensors)  and 
Dmax  =  L/c.  When  the  source  is  in  the  broadside  position  (along  the 
y-axis  in  Figure  1.2)  the  time  delay  is  zero.  The  time  delay  is  known 
"a  priori"  to  be  limited  by  -L/c  <  D  <  L/c,  or  equivalently,  -L  <  aR  <  L. 


Many  techniques  have  been  proposed  for  TDE  including  both  frequency 
domain  and  time  domain  processors.  Perhaps  the  simplest  method  for 
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estimating  the  time  delay  is  to  compute  the  cross-correlation  function 
between  the  signal  and  its  del ayed  version.  The  argument  for  which  the 
cross-correlation  function  attains  its  maximum  value  corresponds  to  the 
time  delay  estimate  [5,  pp.  64-65].  Knapp  and  Carter  have  shown  that  many 
of  the  commonly  used  TDE  techniques  are  related  through  the  Generalized 
Cross  Correlation  (GCC)  approach  [6].  The  class  of  GCC  processors  are 
generally  considered  to  be  frequency  domain  processors  and  include;  the 
Roth  processor  [7],  the  Eckart  filter  [8],  the  smoothed  coherence  transform 
(SCOT)  [9],  the  phase  transform  (PHAT)  [10],  the  Wiener  processor  [11],  and 
the  Hannan-Thomson  or  maximum  likelihood  processor  [12,  6].  In  this  work, 
consideration  is  limited  to  the  GCC  approach  for  TDE  which  is  discussed  in 
some  detail  in  Chapter  2.  However,  it  should  be  noted  that  considerable 
work  has  been  done  on  time  domain  processors  for  TDE.  Many  of  these  time 
domain  techniques  are  based  on  Widrow's  least  mean  square  (LMS)  algorithm 
[13].  Excellent  discussions  of  adaptive  time  delay  estimation  are  given  in 
works  by  Youn  [14],  Youn  and  Ahmed  [15],  Chan,  Riley,  and  Plant  [16,  17], 
Etter  and  Stearns  [18],  and  Feintuch,  Bershad,  and  Reed  [19]. 

1.2  TDE  Performance 

Several  of  the  GCC  processors,  such  as  the  maximum  likelihood  and 
Eckart  estimators,  are  designed  to  be  optimum  with  respect  to  some 
performance  criteria.  Others,  such  as  the  SCOT  and  PHAT,  are  intuitive  or 
"ad  hoc"  techniques,  developed  to  perform  well  for  certain  signal  and  noise 
spectra.  Depending  upon  the  performance  criteria  and  the  signal  and  noise 
spectra,  a  strong  argument  can  be  made  for  the  optimality  of  each  of  the 
TDE  processors.  There  is  an  obvious  need  for  a  performance  standard 
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against  which  the  various  TDE  methods  can  be  evaluated.  The  Cramer-Rao 
lower  bound  (CRLB)  is  commonly  used  as  a  performance  standard  [20-22].  The 
CRLB  yields  a  lower  bound  on  the  variance  of  any  unbiased  time  delay 
estimate  as  a  function  of  several  parameters  (e.g.  the  signal  and  noise 
power  spectra  and  the  observation  time).  Part  of  the  appeal  of  the  CRLB  is 
that  under  certain  assumptions,  there  is  a  theorem  which  states  that  the 
maximum  likelihood  (ML)  estimate  can  be  made  arbitrarily  close  to  the  CRLB 
for  sufficiently  long  observation  times  [23,  pp.  62-72].  However,  the 
theorem  does  not  specify  how  long  the  observation  time  must  be.  Thus, 
while  the  CRLB  sets  a  lower  bound  on  the  variance  of  the  time  delay 
estimate,  actual  performance  can  be  much  worse  for  a  given  signal-to-noise 
ratio  (SNR)  and  observation  time.  This  is  substantiated  by  the  simulation 
results  in  [24,  25]. 

Several  studies  have  been  conducted  to  find  a  bound  tighter  than  the 
CRLB,  which  would  predict  performance  more  accurately.  Chow  and 
Schultheiss  have  investigated  a  simplified  version  of  the  Barankin  bound 
[26].  lanniello  has  developed  a  correlator  performance  estimate  (CPE)  and 
has  shown  via  simulation  that,  for  the  cross-correlation  technique  of  TDE, 
the  CPE  yields  a  more  accurate  estimate  of  performance  than  the  CRLB  [27]. 
Weiss  and  Weinstein  have  proposed  a  modified  version  of  the  Ziv-Zakai  lower 
bound  (ZZLB)  [28-30].  A  comparison  of  the  CPE  and  the  ZZLB  is  reported  by 
lanniello,  Weinstein,  and  Weiss  in  [31].  Further  comparisons  of  the  CPE, 
the  ZZLB  and  the  CRLB  are  presented  in  this  dissertation.  In  particular, 
the  behavior  of  these  performance  estimates  is  considered  as  a  function  of 
the  observation  time  and  SNR,  and  the  implications  of  this  behavior  are 
discussed.  Portions  of  this  work,  comparing  the  CRLB  and  the  CPE,  appear 
in  [32]. 
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1.3  Organization  and  Content 


Before  proceeding  with  the  analysis  of  time  delay  estimator 
performance,  the  mathematical  model  for  the  TOE  problem  is  introduced  in 
Chapter  2.  The  GCC  approach  for  TOE  is  discussed  briefly  and  several 
specific  estimators  from  the  class  of  GCC  processors  are  described.  In 
Chapter  3,  the  shortcomings  of  the  CRLB  are  pointed  out,  in  order  to 
motivate  the  need  for  better  performance  estimates  such  as  the  CPE  and  the 
ZZLB.  Derivations  of  the  CRLB,  the  CPE  and  the  ZZLB  are  outlined  and 
expressions  for  the  variance  of  the  delay  estimate  for  specific  signal  and 
noise  power  spectra  are  given.  Futher  comparisons  of  these  performance 
estimates  are  presented  in  Chapter  4,  including  comparison  with  simulation 
results.  The  threshold  effect  exhibited  by  the  CPE  and  the  ZZLB  is  shown 
to  have  important  implications  for  TDE  processing  methods.  In  particular, 
significant  performance  gains  may  be  obtained  from  increasing  the  coherent 
processing  time  as  opposed  to  increasing  the  incoherent  processing  time. 

To  this  point,  a  stationary  TDE  model  (no  relative  motion  between  source 
and  sensors)  has  been  assumed.  In  Chapter  5,  the  case  of  a  time-varying 
time  delay  is  considered.  The  observations  of  Chapter  4  point  to  the 
necessity  of  pre-processing  the  received  signals  to  compensate  for  the 
relative  motion  before  applying  the  GCC  approach  for  TDE.  One  such 
compensation  technique  is  described  and  preliminary  simulation  results  are 
presented.  Chapter  6  summarizes  the  conclusions  of  this  research  and 
offers  some  suggestions  for  future  work. 
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CHAPTER  2 


THE  GCC  APPROACH  FOR  TIME  DELAY  ESTIMATION 

This  chapter  introduces  the  mathematical  model  and  discusses  the 
related  assumptions  for  the  time  delay  estimation  problem  considered  in 
this  work.  Based  on  this  model,  the  Generalized  Cross  Correlation  (GCC) 
approach  for  estimating  time  delay  is  developed.  Several  specific 
techniques  from  the  class  of  GCC  estimators  are  formulated  and  the 
characteristics  of  these  techniques  are  discussed  briefly. 

2.1  Mathematical  Model 

A  commonly  used  model  for  a  signal  radiated  from  an  acoustic  source 
and  received  in  the  presence  of  additive  noise  at  two  spatially  separated 
sensors  is  given  by  (2-1).  Denoting  the  received  signals  as  r^(t)  and 
r2(t),  let 

r^(t)  =  s(t)  +  n^(t)  (2-la) 

r2(t)  =  as(t+D)  +  n2(t),  0  £  t  £  T  (2-lb) 

where  s(t)  represents  the  source  signal,  n^(t)  and  n2(t)  are  the  additive 
noises,  a  is  a  relative  attenuation  parameter,  and  D  is  the  time  delay 
parameter  to  be  estimated.  The  observation  time,  T,  is  assumed  to  be  much 
larger  than  D.  Additionally,  it  is  assumed  that  s(t),  n^(t),  and  n2(t) 
are  zero  mean,  stationary,  Gaussian  random  processes  and  that  they  are 
mutually  uncorrelated.  Also  note  that  the  model  of  (2-1)  makes  several 
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implicit  assumptions.  In  particular,  the  effects  of  multiple  sources  and 
multi-path  arrivals  are  ignored.  Since  only  the  direct  acoustic  path  is 
considered,  the  problem  is  essentially  planar,  with  the  source  and  the  two 
sensors  defining  the  plane  of  interest,  as  in  Figure  1.1.  The  attenuation 
parameter,  o,  allows  modelling  of  a  relative  attenuation  in  the  signal 
strength  due  to  the  difference  in  path  lengths  between  the  source  and 
sensors.  In  general,  o  will  be  frequency  dependent  rather  than  constant  as 
in  (2-1).  Also  the  time  delay,  D,  is  assumed  to  be  constant  over  the 
observation  interval  which  requires  that  the  source  remain  stationary 
relative  to  the  sensors  during  this  time  interval.  The  problem  of  relative 
motion  between  source  and  sensors  is  discussed  in  Chapter  5.  Finally,  the 
speed  of  sound  in  the  ocean  is  assumed  constant,  although  it  actually 
varies  with  depth  (temperature  and  pressure). 

The  results  of  any  analysis  based  on  the  simplified  model  of  (2-1) 
must  be  applied  cautiously  to  "real  world"  problems,  such  as  the  passive 
sonar  problem  which  motivates  this  research.  Nevertheless  considerable 
insight  into  time  delay  estimator  performance  can  be  derived  from  analysis 
of  this  simple  model.  More  detailed  discussions  of  sonar  signal  processing 
considerations  and  effects  of  the  ocean  medium  may  be  found  in  [33-35]. 

2.2  GCC  Approach 

An  effective  method  of  estimating  the  time  delay,  D,  is  to  compute  the 
cross  correlation  function  between  the  received  signals,  and  r2(t). 
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(2-2) 


-  (t)  =  E[r,(t+T)  r„(t)] 
i  ^ 

where  E[*]  is  the  expectation  operator.  Substituting  (2-1)  into  (2-2) 
yields 


R„  „  (t)  =  a  R„(t-D),  (2-3). 

2 

where  ^^^(t)  is  the  auto-correlation  function  of  the  signal,  s(t).  To 
obtain  (2-3),  the  assumption  that  s(t),  n2^(t),  and  n2(t)  are  mutually 
uncorrelated  has  been  used.  The  cross  correlation  function  can  also  be 
computed  from  the  cross  power  spectrum,  G  (f),  of  the  received  signals, 

''r2 

using  the  Fourier  transform  (FT)  relationship 

oo 

^  r  ^r  r  (2-4a) 

where,  for  the  model  of  (2-1), 

G^^^^(f)  =  a  G^5(f)  .  (2-4b) 

It  can  be  seen  from  either  (2-3)  or  (2-4)  that  R  (x)  attains  its  maximum 

rir2 

value  when  t=D.  Thus,  the  argument,  x,  which  maximizes  the  cross 
correlation  function  is  an  estimate  of  the  time  delay. 

In  practice,  the  cross  correlation  function  must  be  estimated  from 
data  acquired  during  a  finite  observation  time,  T.  To  improve  the  accuracy 
of  the  delay  estimate,  D,  it  is  usually  desirable  to  prefilter  the  noisy 
received  signals,  r2^(t)  and  r2(t),  prior  to  computing  the  cross 
correlation  function.  For  example,  if  the  source  signal  is  known  to  be 
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bandl imited  between  and  f2,  while  the  spectra  of  the  additive  noises 
extend  beyond  these  limits,  it  would  be  reasonable  to  prefilter  the 
received  signals,  passing  only  the  frequency  band  between  f^^  and  f2. 
This  concept  of  prefiltering  motivates  the  GCC  approach  for  time  delay 
estimation  as  depicted  in  Figure  2.1  [6]. 


The  received  signals,  and  r2(t),  are  processed  through 

filters  having  transfer  functions,  H^(f)  and  H2(f),  respectively.  The 
cross  power  spectrum  between  the  filter  outputs,  P]^(t)  and  P2(t),  is 
given  by  [36,  p.  353] 

-  "l'f>  Vz"’’ 


where  *  denotes  the  complex  conjugate.  The  GCC  function,  R  (t),  is 

rir2 

then  obtained  by  taking  the  inverse  Fourier  transform  (IFT)  of  Gp  ^  (f)^ 


r  =f  W(n  „  (f)  df,  (2-6) 


where  W(f)  =  H^(f)  H2(f).  Substituting  (2-4b)  in  (2-6)  yields 

oo 

4  r  df.  (2-7) 

•  '  1'  2  -*00 


0 

To  ensure  that  (x)  has  a  maximum  at  t=D,  W(f)  must  be  a  real 

'^r2 

function,  or  equivalently,  the  prefilters,  H^(f)  and  H2(f),  must  have 
identical  phase  characteristics. 


While  the  realization  of  Figure  2.1  is  useful  conceptually,  the  approach 
suggested  by  (2-6)  is  often  preferred.  That  is,  the  cross  power  spectrum. 
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Figure  2.1  GCC  Approach  for  IDE 


-  (f),  is  multiplied  by  a  real  weighting  function,  W(f),  then  the  IFT  of 

the  result  is  computed  to  obtain  the  GCC  function.  Additional  insight  into 
the  GCC  approach  may  be  gained  by  writing  (2-6)  in  a  slightly  different 
form. 


00 

R®  „  (t)  - /w.(f)  df, 

12  ^ 


(2-8a) 


where 


W^(f )  -  W(f) 


Gr  ,  (f) 

^rz 


a  W(f)  G^^(f),  (2-8b) 


and 


rir2 


G.  .  (f) 


'"l'"2 


e-j2irfD  ^  ^2_8c) 


The  derivative  of  the  phase,  d(f),  with  respect  to  the  frequency  is 
a  measure  of  the  time  delay  and  the  GCC  method  can  be  thought  of  as 
applying  a  weighting  function,  W^(f),  to  the  phase  of  the  cross  power 
spectrum  [37].  Note  that  (2-8)  is  equivalent  to  (2-6). 


A  variety  of  weighting  functions  have  been  proposed  to  optimize 
performance  with  respect  to  certain  criteria.  As  suggested  previously, 
if  a  priori  information  about  the  signal  and  noise  spectra  exists,  this 
information  should  be  incorporated  in  the  weighting  function,  W(f). 
Indeed  many  of  the  commonly  used  weightings  are  functions  of  the  signal 
and  noise  power  spectra.  However,  in  most  cases  of  practical  interest, 
little  or  no  a  priori  information  is  available.  This  necessitates  the 
estimation  of  the  cross  power  spectrum  and  the  weighting  function  from 
the  received  signals.  The  pertinent  spectra  can  be  obtained  using 
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Processor 

W(f)  =  H^(f)  H*(f) 

W^(f)  - 

Gr  r  (f) 
rir2 

W(f) 

Standard  Cross 
Correlation  (SCC) 

'  1 

^  (f) 
12 

Roth 

1 

Gr  r  if) 

G  (f ) 

''ri 


G  (f ) 


Smoothed  Coherence 
Transform  (SCOT) 


2  2 


Wiener  Processor  (WP) 


Cr  .  (f) 

'^r  2 


Cr  r  (f) 


Gr  r  (f) 


Phase  Transform  (PHAT)  G  (f) 


''l''2 


Maximum  Likelihood  (ML) 


''l''2 


[i  -  (f)] 


''l''2 


(f) 

IT 


'"l'"2 


(TT 


^  r 

l~rz — 

’"l’"2 


IT) 


Table  2.1  Weighting  Functions  for  GCC  Processors 


15 


standard  spectral  estimation  techniques.  Substitution  of  these  estimates 
into  (2-6)  or  (2-8a)  enables  an  estimate  of  the  GCC  function, 

^  Q 

-  ('^)»  to  be  computed.  Inherent  in  the  need  to  estimate  the  power 
2 

spectra  is  a  corresponding  degradation  in  performance,  the  degree  of  which 
is  dependent  on  the  merit  of  the  estimates.  This  point  will  be  considered 
again  in  connection  with  the  simulation  results  in  Chapter  4.  In  the 
following  discussion,  specific  techniques  for  time  delay  estimation  are 
examined  and  it  is  shown  how  they  are  related  through  the  GCC  approach. 

The  GCC  weighting  functions  for  these  techniques  are  summarized  in  Table 

2.1. 


2.2.1  Standard  Cross  Correlation 


The  basis  for  the  GCC  approach  is  the  standard  cross  correlation  (SCC) 
method.  As  noted  in  (2-4a),  the  SCC  function  and  the  cross  power  spectrum 
form  a  FT  pair.  The  SCC  method  can  be  thought  of  as  applying  a  uniform 
weighting  to  the  cross  power  spectrum  prior  to  computing  the  IFT.  Setting 


W(f)  =  1  in  (2-6),  or  equivalently,  W^(f)  = 


rir2 


(f) 


in  (2-8a), 


immediately  yields  the  defining  relationship  for  the  SCC  function,  equation 
(2-4a). 


2.2.2  The  Roth  Processor 


One  of  the  first  modifications  to  the  SCC  method  was  proposed  by 
Roth  [7].  The  Roth  processor  applies  the  weighting  function 
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{2-9a) 


to  obtain 

(t)  =  r df.  (2-%) 

If  r^{t)  is  the  input  to  a  linear  filter  which  provides  the  minimum  mean 
square  error  estimate  of  r2{t)  (i.e.  a  Wiener  filter),  then 

R^°th  )  -jg  impulse  response  of  this  filter.  Thus,  the  transfer 
'”l'”2 

function  of  this  optimum  filter  is  simply  G  (f)/G  (f). 

’"r  1 

For  the  signal  model  under  consideration, 

G.  .  (f)  =  G  (f)  +  G  (f)  .  (2-10) 

"r  1 

The  Roth  weighting  function  suppresses  those  frequency  bands  where  the 

noise  power,  G  (f),  is  large  and  the  estimate  G  (f)  is  more  likely 
”l”l  '  r  2 

to  be  in  error.  Further  the  Roth  weighting  tends  to  suppress  those  regions 
where  G^^{f)  is  large,  which  has  a  whitening  effect  on  the  signal.  This 
reduces  the  ambiguity  that  is  introduced  in  the  SCC  function  by  strong 
tonals  in  the  signal  spectrum  [9,  36]. 


2.2.3  The  Smoothed  Coherence  Transform  (SCOT)  ■ 

Unless  it  is  known  that  either  G„  „  (f)  or  G  „  (f)  is  the 

n^n^' 

dominant  noise  process,  there  is  no  reason  to  choose  the  weighting  function 
W(f)  =  1/G^  ^  (f)  over  W(f)  =  1/G^  ^  (f).  The  SCOT  [9]  avoids  this 

riri  r^r^ 
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uncertainty  by  selecting 


W(f)  =  1^ 


(2-11) 


The  SCOT  is  then  given  by 


(2-12a) 


where 


(2-12b) 


In  (2-12b),  Y  (f)  is  the  complex  coherence  between  r,  (t)  and  r,(t). 

rir2  i  ^ 

In  terms  of  the  GCC  realization  of  Figure  2.1,  the  SCOT  can  be 

interpreted  as  employing  two  pre-whitening  filters,  H, (f)  =  l/yo (f) 

i  r.  r. 


and  H2(f)  =  1/v/G^  ^  (f ) ,  followed  by  a  cross  correlator.  Thus  the  SCOT 


2' 2 


exhibits  advantages  similar  to  the  Roth  processor,  suppressing  frequency  bands 
where  the  noise  power  is  large  and  reducing  the  effect  of  strong  tonal s  in  the 
signal,  while  eliminating  the  ambiguity  of  the  Roth  processor  as  to  which 
noise  process  is  dominant.  If  G  (f)  =  G  (f),  the  SCOT  and  the  Roth 

technique  are  identical.  Writing  the  SCOT,  equation  (2-12a),  in  the  form  of 
(2-8a)  yields 


C„  „  (f)  df  (2-13) 


where  C  (f)  is  the  magnitude  squared  coherence  (MSC),  that  is 

rir2 
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2 


(2-14) 


r  = 

r2 


r 


The  SCOT  assigns  greater  weight  to  the  phase  (time  delay)  information  in 
frequency  bands  where  the  coherence  is  high. 


2.2.4  The  Wiener  Processor 


Recently  a  new  GCC  processor,  referred  to  as  the  Wiener  processor  (WP) 
by  the  authors,  has  been  proposed  in  [11].  This  processor  uses  two  Wiener 
prefilters  to  obtain  the  minimum  mean  square  estimates  of  s(t)  from  rj^(t), 
and  of  as(t-D)  from  r2(t)  before  computing  the  cross  correlation.  The 
prefilters  are  given  by 


and 


G  (f ) 
ss'  ' 


""ri 


G  (f) 

1  '"l'"2  j2irfD 

a  GTTTfl" 

1 


(2-15a) 


a%s(^) 


*^2^^^  ■  G  ■ 


'"2'"2 


a  G  (f) 

'^l'^2  j2irfD 

G - (fT"  ® 

r  r  ' 


(2-15b) 


The  corresponding  weighting  function  is  then  the  MSC, 


W(f)  =  H  (f)  H  (f)  =  ^  (f) 

i  r^r^ 


(2-16a) 


and  the  WP  computes  the  GCC  function, 

nWP 


rir2 


^.=7 


.  „  (f)  G„  „  (f)  df. 

1*^2  ’^l’^2 


(2-16b) 


The  WP  suppresses  the  cross  power  spectrum  in  frequency  regions  where  the 
coherence  is  low.  Note  that  while  the  Wiener  prefilters,  Hj^(f)  and  H2(f), 
require  knowledge  of  the  attenuation  parameter,  a,  the  weighting  function 
can  be  estimated  directly  from  the  received  signals. 
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2.2.5  The  Phase  Transform  (PHAT) 


The  PHAT  is  an  "ad  hoc"  technique  which  uses  the  weighting  [10] 


or 


W(f)  =  1/ 


=  1. 


^  r 

rir2 


(2-17a) 


(2-17b) 


Thus  the  PHAT  computes  the  IFT  of  the  phase  function. 


r  1  r*  o 
12 


where  for  the  model  of  (2-1) 

^(f)  =  -2TrfD 

so  that,  ideally 

,PHAT 


R 


'"l'"2 


(t)  =  6(t-D)  . 


(2-18a) 


(2-1 8b) 


(2-1 8c) 


In  practice,  only  an  estimate,  ^(f),  of  the  phase  function  can  be 
PHAT 

obtained,  and  R„  (t)  will  not  be  an  ideal  delta  function.  The  PHAT 

rir2 

assigns  equal  weighting  to  the  phase  estimate  throughout  the  frequency 
band,  independent  of  SNR.  Thus,  the  PHAT  fails  to  suppress  frequency  bands 
where  the  SNR  is  relatively  low  andwhere  the  phase  estimate  is  more  likely 
to  be  in  error.  However,  the  PHAT  effectively  whitens  the  cross  power 
spectrum  and  therefore  eliminates  the  effect  of  strong  tonal s. 


2.2.6  The  Maximum  Likelihood  (ML)  Estimator 


The  final  weighting  to  be  considered  is  that  of  the  ML  estimator  [6]. 
This  technique  is  identical  to  that  proposed  by  Hannan  and  Thomson  [12]  and 
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is  often  referred  to  as  the  HT  processor.  The  ML  designation  will  be  used 
in  this  work.  As  shown  in  Appendix  A,  the  ML  estimator  is  obtained  by 
applying  the  weighting  function 


thus 


W(f)  = 


^r2 


ITT 


Cr  ^  if) 

r--c;7 

*^1^2 


m 


(2-19a) 


^ML 

^1^2 


Cr  r  (f) 
r,r2 


P  (f) 


(2-19b) 


The  ML  estimator,  like  the  SCOT,  weights  the  phase  relative  to  the 

c  pnT 

MSC.  However,  (f)  =  y/C  (f)  approaches  unity  when  the  coherence  is 

0  '  r  2 

large  (note,  0  <  ^  (f)  <  1),  whereas  wJJ^f)  =  ^  (f)/(l  -  ^  (f)) 

^12  ^  12  12 
goes  to  infinity  as  the  MSC  approaches  one.  Compared  to  the  SCOT,  the  ML 

estimator  assigns  much  greater  weight  to  the  phase  in  frequency  bands  where 

the  coherence  is  large.  While  the  SCOT  reduces  the  effect  of  strong 

tonals,  the  ML  estimator  tends  to  emphasize  such  tonals,  where  the 

coherence  will  be  relatively  large.  In  practice,  the  pertinent  power 

spectra  must  be  estimated  and  only  an  estimate  of  (2-19)  can  be  obtained. 

Hence,  only  an  approximate  maximum  likelihood  (AML)  estimator  can  be 

implemented.  .  . 
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CHAPTER  3 


TIME  DELAY  ESTIMATOR  PERFORMANCE  BOUNDS 


The  Cramer-Rao  lower  bound  (CRLB)  is  commonly  used  to  set  a  bound  on 
the  variance  of  the  time  delay  estimate.  Use  of  the  CRLB  to  predict  time 
delay  estimator  performance  is  justified  in  that  the  ML  estimator 
performance  approaches  CRLB  performance  for  sufficiently  long  integration 
times,  that  is,  the  ML  estimator  is  asymptotically  efficient  [23,  p.  71]. 
In  practice,  the  observation  time  is  necessarily  limited,  and  in  many 
cases,  the  observation  time  will  not  be  "sufficiently  long"  to  attain  CRLB 
performance.  This  leads  to  two  important  questions  in  the  TDE  problem. 
First,  under  what  conditions  is  CRLB  performance  attainable,  and  second, 
what  performance  can  be  expected  when  conditions  do  not  allow  for  CRLB 
performance.  In  this  chapter,  the  CRLB,  a  correlator  performance  estimate 
(CPE)  and  a  modified  Ziv-Zakai  lower  bound  (ZZLB)  are  investigated  to  find 
answers  to  these  equations. 

3.1  The  Cramer-Rao  Lower  Bound  (CRLB) 

Consider  any  unbiased  estimate,  D(R^),  of  the  time  delay  D,  where  ^  is 
an  observation  vector  with  conditional  probability  density  function 
p(R^|t).  Then  the  following  inequality  for  the  variance  of  the  delay 
estimate  holds  [23,  p.  66]  or  [38,  p.  316], 


0?  =  E[(D(R)  -  dfy 
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3  In  p( R 1 x) 
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(3-1) 
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provided  that  8p(r|t)/8t  and  3^p(r|t)/3t^  exist  and  are  absolutely 
integrable.  From  equation  (A-8),  (see  Appendix  A),  it  follows  that 
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3  In  p ( R  I  t) 


3t 
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=  E 
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(3-2) 


since  the  quantity  c  in  (A-8c)  is  independent  of  t.  Substituting  (A-10) 
into  (3-2)  and  noting  that  the  term  is  also  independent  of  t  yields 
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( 3-3c ) 


To  obtain  (3-3c)  from  (3-3b),  the  relation  E[R. (f)  R«(f)]  =  T  G  (f) 

i  ^  r^r2 

has  been  used.  Substitution  of  (3-3c)  into  (3-1),  after  taking  the 
indicated  second  partial  deri vativ£. with  respect  to  t ,  results  in  the 
following  expression  for  the  CRLB, 
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The  CRLB  depends  upon  the  signal  and  noise  power  spectra,  in  the  form  of  the 
magnitude  squared  coherence  (MSC),  and  is  inversely  proportional  to  the 
observation  time,  T. 


Now  consider  specific  signal  and  noise  power  spectra, 

0  <  I  f  I  <  B 


G,s(f)  = 


el sewhere 


(3-5a) 


and 


G  (f )  = 
nn'  ' 


No/2  .  o<|f|<B 


0  ,  elsewhere 


(3-5b) 


where  G  (f)  =  G  (f)  =  G  (f).  In  this  work,  the  signal-to-noi se 
nn  n^n^ 

ratio  defined  as 


SNR  i 


0  ,  elsewhere 


(3-6) 


Note  that  this  is  the  SNR  at  the  input  of  a  single  receiver.  For  these 
power  spectra,  the  MSC  is  given  by 

Gss(f)^ 

(f)  =  - 52 - ^  (3.7a) 


(Gss(f)  -G^(f))- 


SNR' 


(SNR  +  1)‘ 


(3-7b) 


Thus  (3-4)  can  be  written  in  terms  of  the  SNR  and  the  CRLB  becomes 
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for  the  power  spectra  of  (3-5). 


(3-8b) 


As  discussed  in  [21],  the  CRLB  exhibits  different  behavior  at  high  and 
low  SNR.  Specifically,  it  can  be  shown  that, 

2  log  ~  ^  ~  ^°9  ^ 

and 

2  log  ~  ^°9  ^  ^°9  ^ 

2 

where  K  =  log(3/8Tr  ).  Thus,  the  high  and  low  SNR  approximations  of 
log  are  linear  functions  of  log  SNR,  with  the  low  SNR  approximation 

having  a  slope  twice  that  of  the  high  SNR  approximation.  Further  it  is 
easy  to  show  that  the  high  and  low  SNR  approximations  intersect  at 
log  SNR  =  -log  2,  that  is,  when  the  input  SNR  at  a  single  receiver  is 
-3  dB,  independent  of  the  observation  time  and  signal  bandwidth. 


(3-9a) 

(3-9b) 


As  discussed  previously,  the  time  delay  is  known  a  priori  to  be 
bounded  by  -L/c  £  £  L/c,  where  L  is  the  sensor  spacing  and  c  is  the 

propagation  speed  of  sound  in  the  ocean.  Most  TDE  methods  inherently 
utilize  this  knowledge  to  limit  the  range  of  the  delay  estimate.  However, 
the  CRLB  does  not  incorporate  this  information  to  limit  the  variance  of  the 
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2 

delay  estimate  and 


oo 


as  SNR  . 


Another  shortcoming  of  greater 


significance  is  that  when  the  observation  time  is  not  sufficiently  long,  - 
the  CRLB  does  not  accurately  predict  time  delay  estimator  performance.  The 
explanation  of  this  deficiency  has  to  do  with  large  estimation  errors, 
which  is  the  subject  of  the  next  section. 

3.2  Large  Estimation  Errors 

The  ML  estimate  is  a  minimum  variance  estimate,  that  is,  it  attains 
the  CRLB  when  the  observation  time  is  long  enough.  Additionally,  when  the 
estimation  error  is  small,  the  variance  of  the  ML  estimate  can  be  expected 
to  be  close  to  the  CRLB.  Therefore,  the  CRLB  is  a  good  estimator  of  TDE 
performance  when  the  observation  time  is  sufficiently  long  to  ensure  that 
the  estimation  error  is  small.  To  gain  a  better  understanding  of  these 
relationships,  a  comparison  of  simulation  results  with  the  CRLB  is  shown  in 
Figure  3.1.  (See  Appendix  B  for  simulation  details.)  The  signal  and  noise 
power  spectra  are  given  by  (3-5)  with  B  =  100  Hz,  relative  to  a  sampling 
frequency  of  2048  Hz.  In  Figure  3.1(a),  T  =  2  seconds,  and  the  simulation 
results  are  in  good  agreement  with  the  CRLB  for  SNR  ^  -3  dB.  Below  -3  dB, 
the  experimental  variance  deviates  sharply  from  the  CRLB.  In  Figure 
3.1(b),  the  observation  time  has  been  increased  to  8  seconds.  Now  the 
experimental  variance  agrees  with  the  CRLB  for  SNR  ^  again 

deviates  from  the  CRLB  below  this  SNR.  Based  on  these  results,  increasing 
the  observation  time  allows  the  estimator  to  maintain  CRLB  performance  for 
lower  SNR,  but  below  some  threshold  SNR,  estimator  performance  begins  to 
degrade  rapidly  relative  to  the  CRLB.  Apparently,  there  is  a  fundamental 
trade-off  between  SNR  and  the  observation  or  coherent  processing  time.  For 
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LOG  ap  log  (7p 


Figure  3.1  Comparison  of  CRLB  and  SCC  Simulation  Results 
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a  given  observation  time,  CRLB,  or  small  estimation  error,  performance  can 
be  obtained  above  the  threshold  SNR.  Below  this  SNR,  large  estimation 
errors  must  be  taken  into  account  and  the  CRLB  is  no  longer  applicable. 

To  understand  the  source  of  such  large  estimation  errors,  consider  the 
cross  correlation  function  for  two  types  of  signals.  First  consider  a 
narrowband  source  such  that  the  bandwidth,  B,  is  small  compared  to  the 

center  frequency,  f^,  (B/f^  <<  1).  A  typical  cross  correlation 

function  for  such  a  source  signal  in  the  absence  of  noise  is  depicted  in 
Figure  3.2(a).  The  maximum  occurs  at  the  true  delay,  D,  but  adjacent  peaks 
are  nearly  as  large.  To  obtain  CRLB  performance,  the  processor  must  be 
able  to  resolve  this  ambiguity  between  peaks,  which  requires  either  a  very 
large  SNR  or  a  very  long  observation  time.  A  similar  effect  is  observed 
for  broadband  signals.  Consider  a  low-pass  signal  as  in  (3-5),  Figure 
3.2(b)  shows  a  typical  cross  correlation  function  for  this  type  of  signal. 
In  this  case,  adjacent  peaks  are  more  widely  separated,  but  at  low  SNR  the 
peak  value  of  the  estimated  cross  correlation  function  may  occur  at  one  of 

these  side  lobes,  yielding  an  estimate  of  the  time  delay  greatly  different 

from  the  true  delay.  These  large  estimation  errors  are  often  referred  to 
as  ambiguities  for  narrowband  signals  and  as  anomalous  estimates  for 
broadband  signals.  In  the  remainder  of  this  chapter,  methods  of  accounting 
for  such  large  estimation  errors  in  the  time  delay  estimate  are 
investigated. 

3.3  The  Barankin  Bound 

As  noted  in  [23,  p.  72],  the  Barankin  bound  is  a  greatest  lower 
bound.  However,  computation  of  the  Barankin  bound  proves  to  be  exceedingly 
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A)  CROSS  CORRELATION  FOR  NARROWBAND  SIGNAL 


B)  CROSS  CORRELATION  FOR  BROADBAND  SIGNAL 


Figure  3.2  Sources  of  Large  Estimation  Errors 
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difficult.  Chow  and  Schultheiss  have  investigated  a  simplified  version  of 
the  Barankin  bound  to  address  the  problem  of  ambiguous  estimates  for  narrow 
band  signals  [26].  While  this  procedure  no  longer  yields  a  greatest  lower 
bound,  the  resulting  lower  bound  has  been  used  with  some  success  to  address 
the  ambiguity  problem.  In  particular,  the  simplified  Barankin  bound 
exhibits  a  distinct  threshold  SNR.  Above  this  SNR,  the  variance  of  the 
time  delay  estimate  is  characterized  by  the  CRLB  and  below  this  SNR  the 
Barankin  bound  exceeds  the  CRLB  by  a  factor  proportional  to  the  square  of 

2  2 

the  ratio  of  the  center  frequency  to  the  bandwidth  (i.e.,  )• 

However,  the  transition  from  the  CRLB  to  the  Barankin  bound  at  the 
threshold  SNR  can  not  be  clearly  specified. 

Also,  like  the  CRLB,  the  Barankin  bound  ignores  the  a  priori 
information  about  the  range  of  possible  delay  values,  and  the  variance  goes 
to  infinity,  as  the  SNR  goes  to  zero.  While  this  approach  succeeds  in 
predicting  the  threshold  phenomena,  it  is  not  the  complete  solution  to  the 
problem. 

3.4  The  Correlator  Performance  Estimate  (CPE) 

The  CPE  was  developed  by  lanniello  [27]  and  provides  an  estimate  of 
the  performance  of  the  SCC  method  of  TDE  (see  Section  2.2.1).  While  this 
approach  does  not  yield  a  lower  bound  on  time  delay  estimator  performance, 
it  does  accurately  estimate  the  performance  of  the  commonly  used  SCC 
method.  The  CPE  computes  the  probability  of  an  anomalous  estimate,  and 
based  on  this  probability,  the  variance  of  the  time  delay  estimate  is 
determined.  The  analysis  presented  here  follows  that  of  lanniello  [27], 
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which  in  turn  is  based  upon  an  approximate  analysis  related  to  pulse 
position  modulation  systems  [39,  pp.  623-640]. 

The  first  step  is  to  define  an  anomalous  estimate  more  precisely. 

Consider  the  SCC  function,  R  (t),  of  equation  (2-3).  Define  the 

rir2 

signal  correlation  time  as 

oo 

Tc4r;^(o)/r,3(t)  d.  .  (3-10) 

—  oo 

The  signal  is  highly  correlated  over  the  correlation  time,  T^,  and  the 
correlation  is  relatively  small  beyond  this  value.  If  the  range  of 
possible  delay  values  is  limited  by  +0^^,  there  are  approximately 

M  =  2D  /T  independent  values  of  the  SCC  function.  Let  R  =  R  (t  ) 

0  c  m 

denote  these  M  independent  values,  and  assume  the  true  delay  is  located 

at  one  of  the  gay  =  D.  An  anomalous  estimate  occurs  when 

|D  -  d|  >  that  is  when  the  delay  estimate  is  farther  than  T^/2 

from  the  true  delay.  The  event  A,  defined  as 

A  =  [R  >  R  ,  for  at  least  one  x  ]  ,  (3-11) 

is  a  reasonable  approximation  to  what  is  meant  by  an  anomalous  estimate 
[38,  p.  629].  Since  the  R^  are  assumed  to  be  independent,  the 
probability  of  A  can  be  formulated  analogously  to  the  probability  of 
error  in  the  communication  of  M  equally  likely  signals,  so  that 
[39,  p.  258] 
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M-1 


dR  (3-12) 
0  '  ' 


where  p(Rq)  is  the  probability  density  function  (pdf)  of  R^,  and 
p(Rj^)  is  the  pdf  for  any  R^^,  which  are  assumed  to  have  identical 
distributions. 

Now  expressions  for  the  pdf  of  R^  and  R^  are  required.  Assume  the 
signal  and  noise  power  spectra  are  as  in  (3-5).  Then  the  autocorrelation 
functions  ^^^(t)  and  ^^^(t)  can  be  written  as 


R33(t)  =  SqB  p(t) 


(3-13a) 


and 


(3-13a) 


where 


sin2TrT  B 
ZittB 


(3-13c) 


is  the  normalized  autocorrelation  function  of  both  the  signal  and  noise. 
Note  that  ^^^(o)  =  S^B  and  Rpp(o)  =  N^B  are  the  variances  of  the 
signal  and  noise,  respectively.  For  a  large  observation  time  -  bandwidth 
product  (BT  >>  1),  the  mean  values  and  variances  of  R^  and  R^  are  given 


by  [40,  p.  183] 


(3-14a) 


4  =■  *  s2]  /K  , 
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(3-14b) 
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and 


cl  s  B^(Sq  +  Nq)^/K  ,  (3-14c) 


where 


K  = 


(3-14d) 


Assuming  that  and  have  Gaussian  distributions,  their  pdf's  are 
completely  specified  by  their  means  and  variances  in  (3-14).  Substitution 
of  the  expressions  for  p(Rq)  and  p(R|^)  into  (3-12)  yields. 
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Making  the  change  of  variable  x  =  R  /a^  and  y  =  »  (3-15) 

0  m 

can  be  written  as 
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and  where  SNR  =  S  /N  .  Evaluating  (3-10)  yields  T  =  1/(2B),  thus,  the 
u  u  ^ 

number  of  independent  values  of  the  SCC  function  is  M  =  D^B.  Equation 
(3-16a)  must  be  evaluated  numerically  to  obtain  the  probability  of  anomaly 
for  a  given  set  of  parameters,  B,  T,  0^  and  SNR,  (see  Appendix  C  for 
computational  details). 

The  variance  of  the  time  delay  estimate  can  be  computed  as  the 
probability  of  anomaly  times  the  variance,  given  an  anomaly,  plus  the 
probability  of  no  anomaly  times  the  variance,  given  no  anomaly.  If  no 
anomaly  occurs,  the  estimation  error  is  small  and  the  CRLB  applies.  When 
an  anomaly  does  occur,  it  has  been  assumed  that  the  anomalous  estimate  is 
equally  likely  to  occur  at  any  of  the  Therefore,  the  variance,  given 

no  anomaly,  is  approximately  that  of  a  random  variable  uniformly 
distributed  in  the  interval  [-0^,  D^],  or  0^/3.  Then,  the  CPE  for 
the  variance  of  the  time  delay  estimate  is  given  by 

4pe  =  PCA]  0^/3  +  (1  -  P[A])  a^pLB  ' 

The  CPE  and  the  CRLB  are  compared  in  Figure  3.3.  The  CPE  is 
characterized  by  three  regions:  1)  at  low  SNR,  P[A]  a  1  and  the  variance 
is  limited  by  prior  information  about  the  range  of  possible  delay  values; 
2)  at  intermediate  SNR,  there  is  a  transition  from  the  prior  information 
limit  to  the  CRLB,  and  3)  at  high  SNR,  P[A]  =  0  and  the  CPE  coincides  with 
the  CRLB.  The  SNR  at  which  the  CPE  begins  to  deviate  from  the  CRLB  is 
referred  to  as  the  threshold  SNR,  (SNR^^^  in  Figure  3.3). 
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Figure  3.3  Comparison  of  CPE  with  CRLB 
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Figure  3.4  Comparison  of  CPE,  CRLB  and  SCC  Simulation  Results 
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The  CPE  is  compared  with  simulation  results  in  Figure  3.4(a)  for 
T  =  2  seconds  and  in  Figure  3.4(b)  for  T  =  8  seconds.  The  theoretical 
variance  is  in  close  agreement  with  the  experimental  variance  and  predicts 
the  threshold  SNR  within  =  1  dB.  Although  the  CPE  is  based  on  an  approximate 
•analysis  and  several  simplifying  assumptions  are  made.  Figure  3.4  shows  that 
the  CPE  is,  nevertheless,  a  very  useful  technique  for  predicting  time  delay 
estimator  performance.  However,  the  CPE  was  developed  specifically  for  the 
see  method  of  TDE  and  is  only  a  performance  estimate,  not  a  performance 
bound.  This  leads  one  to  question  whether  there  may  be  a  time  delay 
estimator  which  can  do  better  than  the  CPE  predicts.  .Therefore,  it  is  still 
of  interest  to  find  a  lower  bound  for  TDE,  which  is  tighter  than  the  CRLB. 
Note  that  if  a  lower  bound  could  be  found,  which  is  close  to  the  CPE,  this 
would  imply  that  the  SCC  method  is  nearly  optimal. 

3.5  The  Ziv-Zakai  Lower  Bound  (ZZLB) 

The  ZZLB  was  derived  by  Ziv  and  Zakai  in  [41],  for  the  signal  parameter 
estimation  problem  in  communications.  Improved  versions  of  this  bound  were 
proposed  independently  by  Chazan,  Zakai  and  Ziv  [42]  and  by  Bellini  and 
Tartar  [43].  Recently,  the  ZZLB  has  been  applied  to  the  time  delay 
estimation  problems  in  works  by  Weiss,  Weinstein  and  lanniello  [28-31].  In 
particular,  the  ZZLB  for  the  variance  of  the  time  delay  estimate  is  derived 
for  the  case  of  narrowband  signals  in  [29]  and  for  the  case  of  low  band 
signals  in  [30].  In  this  section  the  ZZLB  is  derived  for  the  broadband 
signals  of  (3-5),  following  the  procedure  outlined  in  [30]. 
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The  ZZLB  is  based  on  the  probability  of  error  in  deciding  between  two 
hypothesized  values,  say  t  and  t  +  a,  (a  >  0),  of  the  parameter  to  be 
estimated.  Here  the  parameter  of  interest  is  again  the  time  delay,  D.  Let  6 
be  an  estimate  of  the  true  delay,  obtained  from  an  arbitrary  estimation 
scheme.  Now  consider  the  following  suboptimal  decision  procedure: 


if|D-T|<|D-T-A|  ,  then  D  =  x, 

if|D-x|>|D-T-A|  ,  then  D  =  x  +  a.  (3-18) 


That  is,  the  decision  rule  selects  x,  if  the  estimate  is  closer  to  x  than  to 
T  +  A,  and  selects  t  +  a,  if  D  is  closer  to  t  +  a  than  to  x.  Note  that  the 
optimum  decision  rule  would  be  based  on  the  likelihood  ratio  test  as 
discussed  in  [23,  pp.  23-27].  The  probability  of  error  for  the  decision  rule 
of  (3-18)  is  given  by, 

P[x]  P[  deciding  t  +  a|t]  +  P[t+a]  P[deciding  x  |  x  a]  (3-19) 

where  P[x]  denotes  the  probability  of  x  occuring  and  P[deciding  x  +  a|x  ]  is 
the  conditional  probability  of  deciding  that  x  a  is  true,  given  that  x  is 
true.  Assuming  that  x  and  x  a  are  equally  likely  to  occur,  (3-19)  can  be 
written  as 

(1/2)  P[D  -  X  >  a/2 |x  ]  +  (1/2)  P[D  -  X  -  A  <  -  a/2  I  X  +  a]  (3-20) 
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Defining  the  estimation  error  as  e  =  D  -  D,  the  minimum  attainable 

probability  of  error  for  this  binary  decision  problem,  x  +  a),  must 

be  less  than  or  equal  to  the  probability  of  error  for  the  suboptimum  decision 
rule  of  (3-18) ,  thus 

Pg(T,  T  +  a)  <  (1/2)  P[c  >  a/2|t  ]  +  (1/2)  P[c  <  -a/2  I  t  +  a].  (3-21) 


The  range  of  possible  delay  values  is  again  assumed  to  be  bounded  by 
Dq,  so  that  x  and  a  must  satisfy  the  condition  that  t,  t  +  a  e  [-Dq,  Dg], 
or  equivalently, 

-Dq  £  X  £  D^-a,  0  £  a  £  2  D^.  (3-22) 

Integrating  (3-21)  with  respect  to  t  over  the  interval  [-Dg,  Dg-A]  yields 


D  -A 


X  +  a)  dx 


-D 


£  (1/2)  y  [p[c  >  a/2|  T  ]  +  p[c  £  -  a/2  I  T  +  a]]  d- 
-D 


D  -A 


=  (!/.)/ 


P[e  >  a/2|x  ]  dx+ 


-D, 


D 

r° 

(1/2)  /  P[e  £  -a/2|t  ]  dx 

-D  +A 
0 


P[|e|  >  a/2|x  ]  dx 


(3-23) 
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Now  let  F(x)  represent  the  average  of  PCUl  _>  x  |t  ],  where  t  i 
distributed  over  [-D^,  D^],  then 


F(x)  4 


P[|  €  1  >  X  |t  ] 


d  T 


=  P[U  I  ^  x] 

0 


and  from  (3-23) 

-A 

Pg(T,  T+  A)  dx  <  F(a/2)  . 

Multiplying  both  sides  of  (3-24b)  by  and  integrating  with 
A  over  [o,  2D|^]  yields 


2D„  D„-a  2D^ 

^  f  A  y*  Pg(T,  T  +  a)  dx  dA  <  J 


aF(a/2)  dA 


0  -D, 


rO 

=  4/  xF(x)dx  ,x  =  a/2 
0 


=  2x^  F(x) 


0 


uniformly 


(3-24a) 


(3-24b) 

respect  to 


(3-25) 
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Since  F(x)  =  P[  |e  |  _>  x],  it  can  always  be  assumed  that  F(2Dq)  =  0.  Note 
if  F(2D^)  ^  0,  this  implies  that  there  is  a  non-zero  probability  that  the 
estimation  error  U  I  >  20^,  or  D  is  more  than  20^  from  the  true  delay,  D. 
However,  it  is  known  a  priori  that  D  e  [-0^,  Dq].  Thus  if  F(2Dq)  ^  0 
the  estimate  can  be  improved  by  an  obvious  modification  (e.g.,  let  D  =  0  if 
|d1  >  D^).  Also  note  that  the  mean  square  error,  or  for  an  unbiased 
estimator,  the  variance  of  the  time  delay  estimate  is  given  by 


2 

D 


dF(x)  . 


(3-26) 


The  ZZLB  of  [42]  is  obtained  by  substituting  (3-26)  into  (3-25)  and 

2 

solving  for  a^. 


2 
a  ^ 

D 


T  a)  dx  dA 


(3-27) 


In  [43],  Bellini  and  Tartara  note  that  F(x)  is  a  non-increasing  function 
of  X.  Therefore,  a  tighter  lower  bound  for  the  right  hand  side  of  (3-24b)  is 
given  by 


G 


T  +  a)  dx 


<  F(a/2) 


(  3-28) 


where  G[*]  is  a  non-increasing  function  of  a,  obtained  by  filling  in  the 
valleys,  if  any  exist,  of  the  bracketed  function.  This  "valley  filling" 
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process  is  illustrated  in  Figure  3.5.  Proceeding  as  before,  the  modified 
ZZLB  is  given  by 

2  1  /*  ° 

‘'d  -  ^  J  "  ® 

0 

In  general,  it  is  not  possible  to  obtain  a  closed  form  expression  for 
Pg(T,  T+  a).  For  the  IDE  model  of  (2-1),  the  Chernoff  bound  can  be  used  to 
obtain  a  good  approximation  of  Pg(T:,  t  +  a).  The  calculation  of  this 
approximation  is  given  in  Appendix  D.  Under  a  large  observation 
time-bandwidth  product  (BT  >>  1),  it  is  found  that 

Pg(T,  T  +  a)  =  ’  independent  of  t  (3-30a) 

and 

Pg(a)  =  exp[a(a)  +  b(a)]-Q(v/2b(a))  ,  (3-30b) 

where 


D  -a 


/ 


P 

e' 


T  +  a)dT 


da. 


(3-29) 


=  -if  In  [1  +  Y(f, 


a(a)  =  -Jj  In  [1  +  Y(f,  a)]  df  , 

0 

00 

b(a)  =jJ  Y(f,  a)/(1  +  Y(f,  M)  df  , 
0 


(3-30c) 

(3-30d) 


'(f.  A)  = 


Ggg(f)  sin^Trfa 

„  (1^)  ^  „  (1^)J  +  G„  „  (f)  G„  „  (f) 


n2n2 


'^2'^2 


(3-30e) 


and 


Q(x) 


(3-30f) 
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Additionally,  it  is  shown  in  Appendix  D,  that  P0(a)  is  bounded  by 


Pg(A)  >  exp 


df 

0 


Since  Pg(A)  is  independent  of  the  expression  for  the  ZZLB  in 
can  be  simplified  to  yield 


2 

> 

D 


1 

W 

0 


A  G[(2Dq-A)  Pg(A)]dA  . 


0 


Now  an  expression  for  the  ZZLB  is  derived  for  the  specific 
and  noise  spectra  of  (3-5).  Equation  (3-30e)  can  be  written  as 


Y(f ,a) 


SNR^  sin^irfA 

?  'SnR'  +  '1" 


=  SNR'  sin^irfA 


where 


SNR' 


A  SNR^ 

“  2SNR  +  1 


and  where  SNR  is  defined  in  (3-6).  Substituting  (3-33b)  into  (; 
(3-30d)  yields 


i(a)  =  -T  1n(l  +  SNR'  sin^irfA)  df 


and 


b(A)  =  T  J  SNR'  sin^iTfA/(l  +  SNR'  sin^irfA)  df 
0 


(3-31) 

(3-27) 

(3-32) 

signal 

(3-33a) 

(3-33b) 

(3-33c) 

l-30c)  and 

(3-34a) 

(3-34b) 
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Taking  the  derivative  of  a(A)  and  b(A)  with  respect  to  a,  one  obtains 
after  simplifying, 


a'(A)  .-t/  df 

B- 

^5^=b'(A)=T/  df 

*i  6(77 


d 


Ca(A) 


+  b(A)]  =  T 


(3-35a) 


(3-35b) 


(3-35c) 


where 

a(f)  =  irf  SNR'  sin2irfA  (3-35d) 

6(f)  =  1  +  SNR'  sin^irfA  .  (3-35e) 

With  a  little  thought,  it  can  be  seen  that  for  a  c  [0,  1/2B],  b'(A)  >_  0 
with  equality  holding  only  for  a  =  0.  Then,  since  b(A)  0  for  ail  a, 
b(A)  is  a  positive  monotonicaliy  increasing  function  of  A.  Thus,  the 
term  Q(y2b(A))  of  equation  (3-30b)  is  monotonicaliy  decreasing  for 
A  €  [0,  1/2B].  Similarly,  a'(A)  +  b'(A)  <_  0  for  ail  a  c  [0,  1/2B]  with 
equality  holding  only  for  a  =  0.  "From  equation  (D-9),  a(A)  +  b(A)  £  0 
for  ail  A,  so  that  a(A)  +  b(A)  is  a  negative  monotonicaliy  decreasing 
function  for  a  e  [0,  1/2B].  Thus,  the  term  exp[a(A)  +  b(A)]  in  (3-30b) 
is  also  monotonicaliy  decreasing.  Therefore,  for  a  g  [0,  1/2B],  Pg(A) 
is  a  monotonicaliy  decreasing  function,  and  the  probability  of  error  in 
correctly  deciding  between  the  hypothesized  delays,  t  and  t  +  a, 
decreases  as  the  separation.  A,  increases.  This  is  an  intuitively 
pleasing  observation. 
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Next  consider  the  behavior  of  Pg(A)  for  a  >  1/2B.  Using  the 
change  of  variable  Q  =  f a  in  (3-34a)  and  (3-34b)  yields 


a(A) 


ln(l  +  SNR' 


2 

sin 


(3-36a) 


and 


b(A)  = 


;/ 


SNR' 


sin^ir  ni{l  +  SNR' 


2 

sin  irfi) 


dsi. 


(3-36b) 


2 

The  integrands  in  these  equation  are  functions  of  sin  irfi  which  is  periodic 
in  9.  with  a  period  of  1.  Note  also  that  sin  irn  is  an  even  function.  Then 
defining  a^  =  n/2B,  where  n  is  a  positive  integer,  (3-36)  becomes 


BT 


n/j 


sin^irj^)  dn 

"  "  -ni2 


and 


n/2 

b(A  )  =  ^  / 

"  "  -n/2 


SNR'  sin^irf^  /(I  +  SNR'  sin^ir^)  df^. 


Now  (3-37a)  can  be  solved  for  a(Af|_)  as  follows. 


(3-37a) 


(3-37b) 


a(A  )  =  -BT  /  ln(l  +  SNR'  sin^irJ^)  df2 
-172 


=  -BT  J  ln(l  + 


coswj^' )  dQ' ,  9'  =  2^  (3-38) 
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This  integral  can  be  found  in  integral  tables  [44,  p.  461,  eq'n  709]  and 
after  some  algebra 


a(A^)  =  -2BT  ln[(l  +  v/l  +  SNR')/2] 


(3-39) 


Similarly  b(A^)  can  be  written  as 


1/2 


b(A^)  =  2BT  f  SNR'  sin^irn  /(I  +  SNR'  sin^ufi)  dfi.  (3-40) 


Again  resorting  to  the  integral  tables  [44,  p.  432,  eq'n  343]  and  a 
little  algebra,  it  can  be  shown  that 


b(,J  ,  bT  -1  .  (3-41 

"  v/l  +  SNR' 

Substituting  (3-39)  and  (3-41)  into  the  expression  for  Pg(A)  in  (3-30b) 
gives 


exp 


- 

’ 

-BT 

21  n  ^ 

1  +  1  +  SNR'^ 

\  yi  +  SNR'  -  1 

y 

v/l  +  SNR' 

X  Q 


Jl  +  SNR'  - 

v/l  +  SNR'  / 


A 


(3-42) 
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This  observation  leads  to  a  bound  for  the  quantity  G[(2Dq  -  a)  Pg(A)] 
in  (3-32).  The  procedure  followed  to  obtain  the  bound  is  described 
pictorially  in  Figure  3.5.  A  possible  curve  for  Pg(A)  is  sketched  in 
Figure  3.5(a),  and  the  corresponding  curves  for  (2D^  -  a)  Pg(A)  and 
G[(2Dq  -  a)  Pe(^)]  sketched  in  Figure  3.5(b).  Define  the  rectangular 
gate  function  U(a,  a^)  as 

1  >  A  -1  <  A  <  A 

U(A.  A„)  =  n-1  -  -  n 

0  ,  elsewhere  . 

Then,  as  seen  from  Figure  3.5(c), 


G[(2Dq  -  A)  Pg(A)]  >  (2Dq  -  aJ  Pg  U(A,  Aj  (3-43a) 

>  (2Dq  -  1/(2B)  -  a)  Pg.  (3-43b) 

The  expression  on  the  right  hand  side  of  (3-43a)  is  represented  by  the 
decreasing  step  function  in  Figure  3.5(c),  and  the  linear,  dashed  line 
represents  the  final  expression  in  (3-43b). 


This  bound  will  prove  to  be  sufficiently  tight  except  near  a  =  0, 
where  the  value  of  Pg(A)  cannot  be  closely  approximated  by  Pg.  In  this 
region,  consider  the  bound  in  (3-31)  for  Pg(A), 


"  CO 

..  CO 

if  \ 

Pg(A)  >  exp 

-  7/  Y^f,A)  df 
0 

.  Q 

\Zjf  Y(f,A)  dA^^^ 
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- ^ ^ ^ ^ ^  A 

1/2B  2/2B  3/2B  4/2B 

A)  POSSIBLE  REPRESENTATION  OF  Pe(A) 


Figure  3.5  Obtaining  a  Lower  Bound  for  G[(2Dq  -  a)  Pe(A)] 
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Now,  simplifying  the  arguments  of  exp[-]  and  Q(-)  in  (3-31), 


2T  /^Y(f,A)  da  -  2T  fsm'  sinSfa  df 
0  *^0 


SNR'  (wfa)'^  df 


where 


2  *2 
=  n  a 


SNR'  B^T/3  . 


Similarly, 


T 

I 


f  Y^(f,a)  df  >  SNR'^  (irfa)^  df 


-V  L 


\ 


where 

SNR'^  T/10. 

Substituting  (3-44)  and  (3-45)  into  (3-31)  yields 


Pg(a)  >  exp(-v^a^)  Q(na), 


for  small  a  . 


(3-44a) 


(3-44b) 


(3-45a) 


(3-45b) 


(3-46) 
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Then,  from  (3-435)  and  (3-46),  the  function  G[{2Dg  -  a)  Pg{A)]  is 
bounded  by 


^{2D  -  a)  exp{-v‘^A^)  Q{iia)  ,  o  <  a  <  z 


G[{2Dq  -  A)  Pg{A)]  >  /  (20^  -  Pg 


,  z  £  A  <_  2D^  -  2g 


’  ^^0  2B  -  ^  -  ^^o‘ 


(3-47) 


The  value  for  z  is  yet  to  be  defined.  Note  that  (20^  -  1/{2B)  -  a)  <  o 
for  A  >  2D^  -  1/{2B).  Therefore,  the  bound  is  set  to  zero  in  the  third 
region  in  (3-47). 


The  choice  of  z  is  somewhat  arbitrary,  but  compare  the  expressions 
for  Pg{z)  in  equations  (3-31)  and  (3-42).  Select  z  to  satisfy  the 
condition 


n  z  =  \/2b{Aj^) 


v/1  +  SNR'  -  1 

n/1  +  SNR' 


1/2 


(3-48) 


Thus,  the  arguments  for  Q{*)  are  chosen  so  that  Q  is  continuous  between 
regions  1  and  2  in  (3-47).  Substituting  (3-445)  into  (3-48)  and  solving  for 'z 
yields. 


z 


^/3  [!-(!  +  SNR')“^^^1 

IT  B  L  SNR'  J 


(3-49) 
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As  SNR'  >  «>  ,  z  »  0,  and  as  SNR'  »  0  ,  application  of  L'Hospital's  rule 
shows  that 

z  ji/z  vB)  -  2-;5  3  ;  thus,  0  <  z  <  ^  . 

The  expression  for  the  ZZLB  is  obtained  by  substituting  (3-47)  into 
(3-32)  to  give 


(2D„  -  i) 


I  '1  'll 
exp(-v  A  ) 


Q(nA)  dA 


f 


20^-1 /(2B) 


A  (2D^  - 


1 

TT 


a)  Pg  dA 


(3-50) 


The  two  integral  terms  on  the  right  hand  side  of  (3-50)  will  be  considered 
separately.  The  first  integral  can  be  simplified  as  follows 


-  A) 


exp(-v^A^) 


Q(nA)  dA 


=  y  (2nDg  -  y)  exp(-vV/n'^)  Q(y)  dy,  y=nA. 

n  0 

Now  consider  the  term  using  (3-44b)  and  (3-45b),  yields 


(3-51) 


4  4  9 

=  wr  ^ 


(3-52) 


since  BT  »  1  by  assumption.  Since  the  major  contribution  of  the 
integral  in  (3-51)  is  for  small  values  of  y,  the  exponential  term  can  be 
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approximated  by  1.  Also  note  that  over  the  range  of  the  integral  in 
(3-51) 


2nDQ  -  y  <  n  (20^  -  z)  =  ZnD^  . 


(3-53) 


Then  (3-51)  becomes 


I 


4.4, 


A  (ZDg  -  a)  exp(-v  A  )  Q(tia)  dA 


^  / 


y  Q(y)  dy 


(3-54) 


The  second  integral  in  (3-50)  can  be  immediately  evaluated  to  give 


I 


20^-1 /(2B) 


A  (2Dq  -  1/(2B)  -  a)  Pg  dA 


=  P. 


\  °  ■  IT) 


a2-  a' 

T  7" 


20^-1 /(2B) 


Po  C 
e  0 


A^/3] 


2D, 


,  2D^  » 


1 

77 


=  Pe  [40^/3  -  Dg  -  z^/3] 


4Pe  D^/3  ,  D^  »  2 


(3-55) 


Substituting  (3-54)  and  (3-55)  into  (3-50)  yields 

iZ 


J  y  Q(y)dy  2P  Dq/3 

D  n  0 


(3-56) 


Thus,  a  closed  form  expression  for  the  Ziv-Zakai  lower  bound  for  the 
variance  of  the  time  delay  is  given  by  (3-56). 
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One  final  simplification  can  be  made  by  evaluating  the  integral  term 


in  (3-56)  as  follows,  using  integration  by  parts 


J  y  Q(y)  dy  =  y  -  J  y-  dy 

0  0 


=  Q(nz)/2  -  (1/2) f  y^Q'(y)  dy 


Recall  that 


q(y)  =  (2.)-'''  /  da  . 


Using  Leibniz  rule  yields 


Q'(y)  =  -(2w)-1/2  e-y 


Then  the  integral  term  on  the  right  hand  side  of  (3-57)  becomes 
iz  nz 


\  ^y^^'(y)  dy  =  fy  .  y  e  ^  dy 

0  2  yJZii  0 


2 


nz 


nz 

exp(-Ti^z^/2) 


-yl2 


dy 


nz 


2  yJZ-n 


exp(-n^z^/2)  +  -^  ^(nz)  -  ^  , 


(3-57) 


(3-58) 


(3-59) 


(3-60a) 

(3-60b) 


(3-60c) 
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where  (3-60b)  is  obtained  from  (3-60a)  by  integration  by  parts,  and 
(3-60c)  is  obtained  by  re-writing  the  integral  term  in  (3-60b)  in  terms 
of  d(y).  Finally  substituting  (3-60)  into  (3-57)  and  (3-57)  into  (3-56), 
the  final  expression  for  the  ZZLB  is  given  by 


2 

D 


> 


-  1]  (6(nz) 


exp(-ii^z^/2)  +  j 

v2ir 


+  2  Pg  D^/3 


(3-61) 


p 

where  nz  is  defined  in  (3-48),  n  is  given  by  (3-44b),  and  is 
given  by  (3-42). 


As  SNR  >0,  Pg  >  1/2  (see  (3-42))  and  nz  >  0  (see  (3-48)),  and  the 
ZZLB  reduces  to 


4|_8  =  d2/3  .  (3-62) 

This  is  the  same  expression  obtained  from  the  CPE  for  low  SNR,  when 
P[anomaly]  =  1.  At  low  SNR,  the  variance  is  limited  by  the  a  priori  limit 
on  the  range  of  possible  delay  values.  At  high  SNR,  Pg  >  0  and  nZ  > 
and  (3-61)  becomes 


2  _  1 
<^ZZLB  =  77 


(3-63) 


2 

Substituting  for  n  from  (3-44b)  and  expressing  SNR'  in  terms  of  SNR  using 
(3-33c)  yields 

2  311+  2SNR 

=  77  •  77  '  ^ 


2 

"  ‘^CRLB 


(3-64) 
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At  high  SNR,  the  ZZLB  coincides  with  the  CRLB.  Thus,  the  ZZLB  exhibits 
behavior  similar  to  that  of  the  CPE.  The  ZZLB,  CPE  and  CRLB  are  compared  in 
Figure  3.6.  The  ZZLB  is  characterized  by  a  prior  information  limit,  a 
transition  region,  and  the  CRLB  at  low,  intermediate  and  high  SNR, 
respectively.  The  only  difference  between  the  CPE  and  ZZLB  is  that  the  ZZLB 
is  somewhat  smaller  than  the  CPE  in  the  transition  region  and  predicts  a 
threshold  SNR  that  is  somewhat  lower  than  that  predicted  by  the  CPE.  This  is 
not  surprising  since  the  ZZLB  is  a  lower  bound  and  the  CPE  is  only  a 
performance  estimate.  A  comparison  of  simulation  results  for  the  SCC  method 
of  TDE  with  the  ZZLB  is  shown  in  Figure  3.7(a)  for  T  =  2  seconds  and  in  Figure 
3.7(b)  for  T  =  8  seconds.-  Note,  these  are  the  same  simulation  results  as  in 
Figures  3.1  and  3.4.  The  simulation  results  and  the  ZZLB  are  seen  to  be  in 
good  agreement,  with  the  ZZLB  predicting  the  threshold  SNR  extremely  well. 

A  listing  of-  a  computer  program  which  calculates  the  CRLB,  the  CPE  and 
the  ZZLB  is  given  at  the  end  of  Appendix  C.  Also  included  are  listings  of 
subroutines  which  evaluate  the  probability  of  anomaly  for  the  CPE  and  the  Q 
function  required  for  both  the  CPE  and  the  ZZLB. 
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Figure  3.6  Comparison  of  ZZLB  with  CRLB  and  CPE 
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Figure  3.7  Comparison  of  ZZLB,  CRLB  and  SCC  Simulation  Results 
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CHAPTER  4 


THRESHOLD  EFFECT  CONSIDERATIONS  FOR  TDE 

The  CPE  and  the  ZZLB  both  predict  a  threshold  effect  in  the 
performance  of  TDE  methods.  Above  a  threshold  SNR,  time  delay  estimator 

performance  is  characterized  by  the  CRLB,  while  below  the  threshold, 

« 

performance  degrades  rapidly  relative  to  the  CRLB.  Simulation  results  for 
the  see  method  corroborate  these  predictions  by  the  CPE  and  the  ZZLB,  as 
seen  in  Chapter  3.  In  this  chapter,  the  threshold  effect  is  investigated 
in  greater  detail  and  the  implications  of  this  effect  are  considered  as 
related  to  coherent  and  incoherent  processing  techniques  of  TDE.  Also, 
additional  simulation  results  are  presented  to  further  substantiate  the 
theoretical  predictions  and  to  allow  comparison  of  several  different  GCC 
methods. 

4.1  Behavior  of  the  Threshold  SNR 

The  effect  on  the  ZZLB  of  varying  the  observation  time,  while  keeping 
B  and  D^  fixed,  is  shown  in  Figure- 4.1.  The  threshold  SNR  decreases  as 
the  observation,  or  coherent  processing  time,  T,  is  increased.  However, 
the  value  of  the  variance  of  the  time  delay  estimate  (or  equivalently, 

log  Op,  as  in  Figure  4.1)  at  which  the  ZZLB  begins  to  deviate  from  the 
CRLB  remains  essentially  constant,  independent  of  the  observation  time.  To 
attain  CRLB  performance  at  SNR  below  the  threshold  SNR  for  a  given 
observation  time,  requires  that  the  observation  time  be  increased 
sufficiently  to  lower  the  threshold  SNR  below  the  desired  operating  SNR. 
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Figure  4.1  Effect  of  Varying  Observation  Time  on  ZZLB 
(B  =  100  Hz,  Dg  =  1/16  second,  T  =  2,  8,  32,  128,  512  seconds) 
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Consider  the  CRLB  and  the  ZZLB  as  functions  of  the  SNR,  with  all  other 
parameters  held  constant.  Below  the  threshold  SNR,  the  ZZLB  quickly  becomes 
an  order  of  magnitude  larger  than  the  CRLB.  Thus,  as  suggested  in  [31],  the 
threshold  SNR  can  be  approximated  by  the  SNR  which  satisfies  the  condition 

All?,  =  ^“^Llb  ’ 


2  2 

where  o^ZLB  “^CRLB  (3-61)  and  (3-64),  respectively. 

The  particular  SNR  for  which  (4-1)  holds  is  given  by 


(4-2a) 


where 

F(x)  =  x2q(x) 


(4-2b) 


and  recall  that 


Q(x)  = 


(4-2c) 


In  (4-2a),  F“^(‘)  denotes  the  larger  of  the  two  solutions  for  the 
inverse  of  F(x).  Also,  recall  the. two  definitions  of  signal-to-noise  ratio 
used  in  this  work. 


SNR  = 


G3s(f) 


0  <  |f|  <  B 


(4-3a) 


and 


SNR' 


SNR^ 

■2'SnR  ■+  '1 


(4-3b) 
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Note  that  (4-2a)  is  expressed  in  terms  of  SNR'.  For  SNR  <<  1, 

p 

SNR'  =  SNR  and  for  large  BT  product,  the  threshold  SNR  becomes  very 
small  as  seen  in  Figure  4.1.  Thus,  for  large  BT,  (4-2a)  can  be  written 
in  terms  of  SNR  as 


Given  that  B  and  are  fixed  and  that  BT  is  large,  SNR^^  is  approximately 
inversely  proportional  to  y/f.  Moreover,  F“^(*)  is  a  relatively  weak 
function  of  its  argument  (i.e.,  an  order  of  magnitude  change  in  the  argument 
produces  a  comparatively  small  change  in  F~^(-).  Therefore,  to  a  good 
approximation,  SNR^^  is  inversely  proportional  to  v/^,  for  large  BT. 
Analytical  considerations  leading  to  (4-2a),  given  (4-1),  are  given  in 
Appendix  E. 


The  behavior  of  the  CPE  due  to  increasing  observation  time  is 
essentially  the  same  as  that  of  the  ZZLB  as  illustrated  in  Figure  4.2. 
Nuttall  has  shown  in  [45]  that  the  threshold  SNR  for  the  CPE  can  be 
approximated  by 


where 


(4-5a) 


(4-5b) 


and  where  Q~^(’)  denotes  the  inverse  of  Q(*)-  For  large  BT,  x  is  small  and 
(4-5)  becomes 


SNRth  = 


1 

v5t 


(4-6) 
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Figure  4.2  Effect  of  Varying  Observation  Time  on  CPE 
(B  =  100  Hz,  Dq  =  1/16  second,  T  =  2,  8,  32,  128,  512  seconds) 
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Note  the  similarity  between  (4-4)  and  (4-6).  The  function  Q“^(*)  is 
also  a  relatively  weak  function  of  its  argument  and  the  threshold  SNR  is 
again  approximately  inversely  proportional  to  However,  to  obtain 

SNRth  from  (4-5)  or  (4-6),  the  probability  of  anomaly,  P[A],  for  which 
the  CPE  begins  to  deviate  significantly  from  the  CRLB  must  be  known.  A 
value  of  P[A]  =  5  x  10“®  has  been  found  empirically  to  give  good 
agreement  between  the  threshold  SNR  computed  via  (4-5)  and  that  observed 
in  Figure  4.2. 

Values  of  the  threshold  SNR  for  the  ZZLB,  equation  (4-2),  and  the 
CPE,  equation  (4-5),  for  various  observation  times  are  given  in  Table 
4.1.  Note  that  these  values  are  in  terms  of  SNR  in  dB.  Equation  (4-4b) 
can  be  used  to  obtain  the  corresponding  values  in  terms  of  SNR'. 


CPE 

ZZLB 

T  (seconds) 

SNRth(dB) 

2 

-2.3 

-4.3 

8 

-6.6 

-8.0 

32 

-10.1 

-11.3 

128 

-13.3 

-14.6 

512 

16.4 

-17.7 

Table  4.1 

Calculated  Threshold  SNR 

Values  for  CPE  and  ZZLB 

(B  =  100  Hz 

,  Dq  =  1/16  sec.,  P[A]  = 

5  X  10“®  for  CPE) 
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4.2  Implications  of  the  Threshold  Effect 


The  threshold  effect  has  particular  significance  relative  to  the 
issue  of  coherent  versus  incoherent  signal  processing  for  TDE.  First, 
observe  that  the  parameter,  T,  referred  to  as  the  observation  time  in  the 
expressions  for  the  CRLB,  the  CPE,  and  the  ZZLB,  is  actually  the  coherent 
processing  time.  If  the  observation  time  is  T  seconds,  the  performance 
predictions  of  Chapter  3  assume  that  the  time  delay  estimate  is  obtained  by 
coherently  processing  all  T  seconds  of  data  coherently,  resulting  in  a 
single  estimate  of  the  time  delay. 

An  alternative  method,  which  is  attractive  in  many  TDE  applications, 
divides  the  T  seconds  of  data  into  N  sections,  each  of  length  T/N  seconds. 
The  N  sections  are  processed  individually  to  obtain  N  time  delay  estimates, 
one  estimate  for  each  T/N  second  data  section.  The  N  delay  estimates  are 
then  averaged  to  obtain  a  final  time  delay  estimate  after  T  seconds.  This 
alternative  method  is  referred  to  as  incoherent  processing.  For  the 
incoherent  processor,  the  coherent  processing  time  is  reduced  to  T/N 
seconds  for  each  of  the  N  sections,  although  the  total  observation  time  is 
still  T  seconds. 

Above  the  threshold  SNR,  time  delay  estimator  performance  is 

2 

characterized  by  the  CRLB,  and  is  inversely  proportional  to  the 

coherent  processing  time.  Now  consider  the  performance  of  the  incoherent 
processor.  For  SNR  greater  than  SNRj^j^,  the  threshold  SNR  for  a  coherent 
processing  time  of  T/N,  the  variance  of  the  time  delay  estimates  obtained 
from  the  T/N  second  data  sections  is  proportional  to  N/T.  The  incoherent 
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processor  then  averages  the  N  time  delay  estimates,  which  reduces  the 
variance  by  a  factor  of  N,  assuming  the  N  estimates  are  independent  [36, 
p.  245].  Thus,  for  SNR  >  SNRj^j^,  the  variance  of  the  time  delay  estimate 
for  the  incoherent  processor  is  proportional  to  (N/T)*(l/N)  =  1/T.  This  is 
the  same  performance  attained  by  coherently  processing  all  T  seconds  of 
data  to  obtain  a  single  time  delay  estimate.  However,  for  SNR  <  SNRj^j^, 
the  performance  of  the  incoherent  processor  degrades  rapidly,  while  the 
coherent  processor  maintains  CRLB  performance  over  the  range 
SNRj  <  SNR  <  SNRj^i^,  where  SNRj  represents  the  threshold  SNR  for  a 
coherent  processing  time  of  T  seconds.  For  this  range  of  SNR,  the  coherent 
processor  exhibits  significant  performance  gains  compared  to  the  incoherent 
processor. 

This  effect  is  illustrated  in  Figure  4.3  for  the  case  of  T  =  8 
seconds  and  N  =  4  sections  (B  =  100  Hz  and  =  1/16  second).  The 
coherent  processor  operates  on  data  sections  8  seconds  long,  -that  is,  the 
coherent  processing  time  is  8  seconds.  The  incoherent  method  processes  the 
data  in  4  sections,  each  2  seconds  long,  so  that  the  total  observation  time 
is  again  8  seconds.  The  predicted  performances  for  the  coherent  and 
incoherent  techniques  are  shown  by  the  curves  labelled  ZZLB  (T  =  8, 
coherent)  and  ZZLB  (T  =  8,  incoherent),  respectively.  The  third  curve, 
labelled  ZZLB  (T  =  2,  coherent),  is  the  ZZLB  for  a  coherent  processing  time 
of  2  seconds.  The  variance  of  the  incoherent  method  is  reduced  by  a  factor 
of  4  relative  to  the  performance  of  a  coherent  processor  for  T  =  2,  due  to 
the  averaging  of  the  4  estimates.  This  appears  as  a  constant  offset 
between  the  curves  ZZLB  (T  =  2,  coherent)  and  ZZLB  (T  =  8,  incoherent)  on 
the  log  oq  vs  SNR  (dB)  plot  of  Figure  4.3.  For  SNR  >  SNR2,  the 
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Figure  4,3  Comparison  of  Coherent  and  Incoherent  Processor  Performance 
(T  =  8  seconds,  N  =  4  sections) 
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performance  of  the  coherent  and  incoherent  methods  is  identical  and 
coincides  with  the  CRLB  fotVf  =  8.  '  For  SNR  <  SNR2,  the  variance  of  the 
incoherent  processor  increases  rapidly,  while  the  coherent  processor 
continues  to  attain  CRLB  performance  until  the  SNR  falls  below  SNRg.  In 
this  example,  the  incoherent  processor  begins  to  deviate  from  the  CRLB  at 
SNR2  =  -2  dB,  while  the  coherent  processor  maintains  CRLB  performance  for 
SNR  approximately  4  dB  lower  than  SNR2  (i.e.,  for  SNR  >  SNRg  =  -6  dB). 

In  general,  for  large  BT  the  threshold  SNR  is  approximately 
proportional  to  l/^T,  as  noted  previously.  Thus,  quadrupling  the  coherent 
processing  time  lowers  the  threshold  SNR  by  a  factor  of  2,  or  by  3  dB. 
Similarly  a  factor  of  64  increase  in  the  coherent  processing  time  will 
lower  the  threshold  SNR  by  approximately  9  dB,  see  Table  4.1  and 
Figure  4.1.  Therefore,  coherent  processing  provides  much  better  estimates 
of  the  time  delay  relative  to  incoherent  processing  for  SNR  below  the 
threshold  for  the  incoherent  processor,  SNR-j-yj^,  but  above  the  threshold 
for  the  coherent  processor,  SNRj.  For  the  signals  considered  here  and 
for  large  BT, 

SHRj  =  SNR^^i^  /n/N  (4-7a) 

or,  in  dB, 

SNR^  (dB)  =  SNR^^I^  (dB)  -  5  log  N.  (4-7b) 

For  example,  if  a  100  second  data  segment  were  available,  and  the 
incoherent  processor  broke  the  data  into  100  sections,  each  1  second  long, 
the  incoherent  processor  would  suffer  a  ~  10  dB  loss  in  the  threshold  SNR 
relative  to  the  coherent  processor.  Performance  gains  can  also  be  realized 
by  increasing  the  signal  bandwidth  for  a  fixed  observation  time,  as  seen 
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from  equations  (4-4)  and  (4-6).  In  many  applications,  and  in  the  passive 
sonar  problem  in  particular,  the  bandwidth  of  the  source  cannot  be  altered 
to  improve  performance.  However,  these  results  indicate  that  the  processor 
should  be  designed  to  operate  over  the  full  bandwidth  of  the  signal,  when 
the  signal  and  noise  spectra  are  as  defined  in  (3-5).  Finally,  note  that 
to  ensure  the  effective  independence  of  the  N  estimates  averaged  by  the 
incoherent  processor,  the  length  of  each  data  section  must  be  much  greater 
than  the  delay  value  (T/N  >>  D).  If  the  N  estimates  are  not  independent, 
the  incoherent  processor  suffers  a  further  degradation  in  performance 
relative  to  the  coherent  processor. 

The  performance  gains,  obtained  by  increasing  the  coherent  processing 
time,  are  achieved  at  the  expense  of  increased  complexity  in  the  coherent 
processing  algorithm.  As  the  coherent  processing  time  is  increased,  it  is 
usually  necessary  to  compensate  for  the  effect  of  time  variations  of  the 
time  delay  (e.g.,  due  to  source  motion).  This  requires  a  substantially 
more  complicated  processor  than  the  standard  GCC  approach  for  the  case  of  a 
fixed  time  delay  as  described  in  Chapter  2.  These  concerns  will  be 
discussed  in  Chapter  5.  Before  proceeding  with  this  topic,  some  additional 
simulation  results  pertaining  to  the  stationary  time  delay  case  are 
presented  to  further  corroborate  the  performance  predictions  of  the  ZZLB. 

4.3  Simulation  Results 

Simulation  results  for  the  SCC,  the  Wiener  processor,  the  SCOT,  the 
PHAT,  and  the  AML  techniques  for  TDE  are  compared  with  the  ZZLB  and  CRLB  in 
Figures  4.4  -  4.8.  Results  are  shown  for  T  =  2  seconds  and  T  =  8  seconds. 
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with  B  =  100  Hz  and  =  1/16  second  for  a  sampling  frequency  of 
2048  Hz.  Complete  simulation  details  are  given  in  Appendix  B.  For  the 
signal  and  noise  spectra  used  in  this  simulation,  see  (3-5),  similar 
results  are  obtained  for  each  of  the  processors,  and  close  agreement 
between  the’  ZZLB  and  experimental  results  is  observed.  However,  note  that 
above  threshold  SNR,  the  experimental  variances  for  the  SCOT,  the  PHAT,  and 
the  AML  techniques  actually  fall  below  the  CRLB. 

This  phenomena  can  be  explained  in  terms  of  an  effective  bandwidth, 
Bg,  as  follows.  The  CRLB  and  ZZLB  were  computed  for  B  =  100  Hz.  As 
discussed  in  Appendix  B,  the  signal  and  noise  sequences  were  obtained  by 
processing  white  noise  sequences  through  a  low  pass  Butterworth  filter, 
which  had  a  cut-off  frequency  of  100  Hz.  Although  a  high  order  filter  was 
used,  the  signal  and  noise  spectra  are  not  ideal  and  contain  some  power  at 
frequencies  above  100  Hz.  Since  the  signal  and  noise  spectra  have  the  same 
shape,  the  SNR,  or  equivalently,  the  MSC  will  remain  approximately  constant 
in  some  small  interval  beyond  B  Hz,  say  B  £  f  £  B^,  until  the  signal 
power  is  diminished  sufficiently  so  that  the  SNR  is  effectively  zero. 

The  phase  weighting  function^,  W^(f),  for  the  five  processors 
considered  here  are  shown  in  Figure  4.9.  While  the  SCC  and  the  MSC 
weighting  functions  go  to  zero  above  100  Hz,  the  SCOT,  the  PHAT  and  the  AML 
weighting  functions  do  not  fall  to  zero  until  f  130  Hz.  The  latter  three 
methods  have  an  effective  bandwidth  of  Bg  =  130  Hz  compared  to 
B  =  100  Hz.  Since  the  SNR  is  approximately  constant,  the  phase  (time 
delay)  information  in  this  extended  bandwidth  is  of  the  same  quality  as  in 
the  nominal  bandwidth.  From  (3-9),  it  is  easily  shown  that  increasing  the 
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Figure  4.4  Comparison  of  SCC  Results  and  ZZLB 
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Figure  4.5  Comparison  of  Wiener  Processor  Results  and  ZZLB 
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Figure  4.6  Comparison  of  SCOT  Results  and  ZZLB 


72 


LOG  an  LOGa 


Figure  4.7  Comparison  of  PHAT  Results  and  ZZLB 
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Figure  4.8  Comparison  of  AML  Results  and  ZZLB 


74 


bandwidth  by  a  factor  of  1.3  has  the  effect  of  subtracting  .17  from  log 
Ocrlb*  Reducing  the  bound  by  this  amount  in  Figures  4.6  -  4.8, 
reconciles  the  theoretical  and  experimental  results  above  the  threshold 
SNR.  Also  note  that  the  threshold  SNR  will  be  slightly  reduced  due  to  the 
increase  in  the  effective  bandwidth;  however,  below  the  threshold  SNR,  this 
increase  in  bandwidth  will  have  little  effect  on  the  ZZLB.  This  phenomena 
is  an  artifact  of  the  technique  used  to  generate  the  signal  and  noise 
sequences  and  does  not  have  great  practical  significance.  However,  it  does 
illustrate  the  whitening  effect  of  the  SCOT  and  the  PHAT  techniques,  and 
that  the  AML  method  weights  the  phase  according  to  the  MSC,  while  the  SCC 
and  WP  methods  weight  the  phase  relative  to  the  magnitude  of  the  cross 
power  spectrum  (i.e.,  the  signal  power),  refer  to  Table  2.1. 
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CHAPTER  5 


ESTIMATION  OF  A  TIME-VARYING  TIME  DELAY 


It  has  been  shown  that  significant  performance  gains  are  realized  at 
low  SNR  by  increasing  the  coherent  processing  time.  For  a  constant  or 
fixed  time  delay,  increasing  the  coherent  processing  time  requires  little 
or  no  modification  of  the  GCC  processor  structure.  However,  when  the  time 
delay  varies  with  time,  the  received  signals  must  be  pre-processed  to 
compensate  for  the  changing  time  delay  to  avoid  degradation  in  performance. 
The  degree  of  degradation  is  dependent  upon  the  amount  of  variation  in  the 
time  delay,  but  even  small  rates  of  change  in  the  delay  can  cause  a 
significant  decrease  in  performance  [46,  47].  In  this  chapter,  a  simple 
model  for  a  time-varying  time  delay  is  considered.  Based  on  this  model, 
the  required  changes  in  the  GCC  structure  to  implement  an  ML  processor  are 
discussed.  A  simplified  compensation  scheme  is  introduced  and  preliminary 
simulation  results  are  presented  to  demonstrate  the  effectiveness  of  the 
compensation  technique. 

5.1  Time-Varying  Time  Delay  ModeJ . 

Consider  the  following  modification  of  the  model  of  (2-1), 


r^(t)  =  s' (6^t)  +  n^(t) 


(5-la) 


r2(t)  =  os'(62t+D))  +  n2(t),  0£t£T 


(5-lb) 


where  rj^(t)  and  r2(t)  represent  the  signals  received  at  two  spatially 
separated  sensors,  s'(t)  is  the  source  signal,  nj^(t)  and  n2(t)  are 
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additive  noises,  a  is  an  attenuation  parameter  and  0  is  the  time  delay  of 
interest.  The  parameters  and  B2  account  for  the  time  compressions 
in  the  source  signal  at  the  two  sensors  due  to  relative  motion  between  the 
source  and  sensors.  If  the  source  is  moving  toward  sensor  1  with  a 
relative  speed  of  Vj^,  then  Bj^  =  1  +  Vj^/c,  where  c  is  the  propagation 
speed  of  the  signal.  For  an  acoustic  signal  in  the  ocean,  c»Vj^,  and 
~  assumed  that  Bj^  and  B2  remain  constant  over  the 

observation  interval  T.  Now  let  s(t)  =  s'(B]^t)  then  (5-1)  can  be  written 
as 

rj^(t)  =  s(t)  +  n^{t)  (5-2a) 

r^{t)  =  as(B(t+0))  +  n^{t)  (5-2b) 

where  B  =  &2l^i  is  the  relative  time  compression.  There  are  now  two 
parameters  to  be  estimated,  the  time  delay,  D,  and  the  relative  time 
compression,  b.  It  is  again  assumed  that  s(t),  nj^(t),  and  n2(t)  are 
sample  functions  from  uncorrelated,  zero  mean,  Gaussian  random  processes. 

Note  that  (5-2b)  can  be  written  as 

r2(t)  =  s(t  +  d(t))  +  n2(t)  (5-3a) 

where 

d(t)  =  d't  +  dg  (5-3b) 

and 

d'  =  (B-1)  ,  =  e  0  .  (5-3c) 

In  (5-3),  d'  represents  the  delay  rate  and  d^  represents  the  initial 
delay.  In  general,  the  time-varying  time  delay  can  be  expressed  as 
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(5-4) 


d(t)  =  cJq  +  d't  +  d"t^  +  d‘"t^  +  ..  , 

where  d^,  d',  d"  ...  are  constants.  The  model  considered  here  assumes 
that  higher  order  terms,  such  as  the  delay  acceleration  (d"),  are 
negligible.  This  model  is  adequate  for  many  applications,  but  is  not 
suitable  for  situations  where  the  higher  order  terms  cannot  be  ignored. 

5.2  The  Maximum  Likelihood  Estimate 


The  cross  correlation  function  between  r^(t)  and  r2(t)  of  (5-1) 
is  given  by 

R^l^^(t+T,  t)  =  aR^^((l  -  B)t  +  T  -  bD)  (5-5a) 


=  aRss(T  - 


d't 


-  V 


(5-5b) 


While  r^(t)  and  r2(t)  are  individually  wide-sense  stationary,  they  .ai'e  no 
longer  jointly  wide-sense  stationary,  (i.e.,  R„  „  (t+T,t)  ^  R-  (t)  for 

te[0,  T]).  The  processes  are  only  jointly  stationary  when  8=1,  that  is 
B^  =  B2*  However,  the  maximum  likelihood  estimator  for  the  dynamic  case 
of  (5-2)  can  be  obtained  following  a  procedure  similar  to  that  used  for  the 
stationary  case.  The  derivation  of  the  ML  estimator  for  b  and  D  has  been 
carried  out  by  Knapp  and  Carter  in  [46]  and  is  outlined  in  Appendix  F  for 
convenience.  For  cases  of  practical  interest,  the  signal  and  noise  spectra 
must  be  estimated  and  approximate  ML  estimates  of  8  and  D  are  obtained  by 
maximizing  the  function  J  with  respect  to  b  and  x,  where 
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oo 


J(b.T)  =  J R^(f)  R*(bf  W(f)df, 


(5-6a) 


—  00 


W(f)  = 


(5-6b) 


G  (f)  (1  -  C 


(f))  , 


and 


>"3(t)  =  r2(t/b)/*v/b 


(5-6c) 


The  values  of  b  and  t  which  maximize  J  are  the  ML  estimates  of  e  and  D, 
respectively.  In  (5-6a),  R^(f)  denotes  the  Fourier  Transform  of  r^(t), 
and  in  (5-6c),  r^(t)  represents  a  time  scaled  version  of  r2(t),  where 
the  time  scaling  counteracts  the  effects  of  the  relative  time  compression. 
Note  that 


R3(f )  =  sTh  R2(bf ) 


(5-7a) 


and  (5-4a)  can  be  written  as 


00 


J(b,T)  =  r R^(f)  R*(f)  W(f)df  . 


(5-7b) 


—  00 


The  ML  estimator  can  be  implemented  as  illustrated  in  Figure  5.1.  A 
value  of  b  is  selected,  r2(t)  is  time  scaled  to  obtain  r3(t)  corresponding 
to  b,  the  FT  of  r^(t)  and  r3(t)  is  computed  and  W(f)  is  estimated,  then 
J(b,T)  in  (5-7b)  is  computed.  The  value  of  x  which  maximizes  J  is  the  best 
estimate  of  D  for  the  b  selected.  This  process  is  repeated  varying  the  value 
of  b  to  generate  an  ambiguity  surface  for  J(b,T)  in  b  and  x.  The  values  of  b 
and  X  corresponding  to  the  global  peak  of  the  ambiguity  surface  are  the  ML 
estimates  of  B  and  0.  In  practice,  knowledge  of  the  physical  problem  is  used 
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Figure  5.1  The  ML  Estimator  for  b,  0 
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to  determine  a  priori  limits  for  s,  and  the  number  of  compensations  (trial 
b  values)  is  chosen  to  provide  the  required  resolution  for  s  over  this 
range.  The  ML  estimator  is  computationally  intensive,  requiring  r2(t) 
and  W(f)  to  be  recomputed  for  each  trial  value  of  b.  Then,  the  function 
J(b,T)  must  be  computed  for  each  b  value  and  this  set  of  functions  must  be 
searched  to  find  the  global  peak.  Additionally  it  is  not  immediately  clear 
how  to  obtain  r2(t)  from  r2(t).  Conceptually,  if  r2(t)  is  recorded 
on  analog  tape,  r2(t)  can  be  obtained  using  a  variable  sampling  rate. 
However,  very  fine  adjustments  of  the  nominal  sampling  rate  are  necessary 
and  there  are  obvious  synchronization  problems.  In  many  cases,  only  a 
sampled  version,  r2(k),  of  r2(t)  will  be  available,  and  a  high  resolution 
interpolation  process  will  be  required  to  obtain  r2(k)  =  r2(k/b) 
from  r2(k). 

5.3  A  Simplified  Compensation  Scheme 

As  discussed  in  Appendix  B,  a  common  technique  for  implementing 
the  GCC  processor  is  to  segment  the  T  second  data  record  into  M 
subintervals.  It  is  often  advantageous  to  use  overlapping  segments 
[48];  however,  disjoint  segments  are  assumed  here  to  simplify  the 
following  discussion.  The  fast  Fourier  transform  of  each  segment  is 
computed  and  the  pertinent  power  spectra  are  estimated  for  each 
segment.  These  M  spectral  estimates  are  then  averaged  in  such  a  way 
that  the  signal  components  will  sum  coherently  to  obtain  a  final  y 
estimate  of  the  cross  power  spectrum  and  the  GCC  weighting  function. 

The  inverse  FFT  of  the  product  of  these  estimates  yields  the  estimate  of 
the  GCC  function. 
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For  the  stationary  time  delay  case  (b  =  1),  the  M  spectral  estimates 
can  simply  be  summed,  and  then  divided  by  M  to  obtain  the  final  spectral 
estimate.  When  b  ^  1,  the  received  signals  must  be  pre-processed  to 
compensate  for  the  relative  time  compression,  or  signal  coherence  will  not 
be  preserved  when  the  spectral  estimates  are  summed.  The  compensation 
method  described  here  compensates  for  the  relative  time  compression  in  the 
frequency  domain  and  avoids  the  need  to  re-process  r2(t)  to  obtain  its 
time  scaled  version,  r2(t),  for  each  hypothesized  value  of  B.  This 
greatly  reduces  the  computational  load. 

Consider  the  expression  for  r2(t)  given  in  (5-3) 

r2(t)  =  s(t  +  d't  +  dg)  +  n2(t),  0<t<T.  (5-8) 

Now  partition  T  into  M  disjoint  segments  of  length  T/M.  Define  the  kth 
interval  as 


T 


k 


(k-l)T/M  <  t  <  kT/M 


and  denote  the  midpoint  of  this  interval  as 


(5-9a) 


t,^  =  (k-l/2)T/M.  (5-9b) 

Then  the  time  delay  at  t  =  tj^  is  given  by  d,^  +  dg,  where  d|^  =  d't|^. 

The  key  step  in  the  development  of  this  compensation  technique  is  to  assume 
that  the  time-varying  time  delay  in  the  interval  Tj^  can  be  approximated  by 
the  value  of  the  delay  at  the  midpoint  of  the  interval,  that  is 
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d't  +  +  (Jq  ,  t  G  T^.  (5-10) 

With  this  approximation,  the  signal  model  becomes 


r^(t)  =  s(t)  +  n^(t) 

(5-lla) 

r2(t)  =  a  s(t  +  d,^  +  d^)  +  n2(t) 

(5-llb) 

for  t  G  T|^,  k=l,  M.  Thus,  r^(t)  and  r2(t)  are  again  jointly 

stationary  over  each  interval  Tj^,  and  the  short  term  cross  correlation 
function  for  this  interval  is  given  by 

Denoting  the  Fourier  transform  for  r.(t),  t  g  T^,  as  R-(f),  where 

I  K  1 

i=l,  2;  the  cross  power  spectrum  for  the  kth  data  segment  can  be  estimated  as 
gI!  „  (f)  =  Ri(f)  R2(f)* 

=  a  G  (f)  .  (5-13) 

s  s-'  ■ 

The  cross  power  spectrum  for  the  kth  interval  undergoes  a  phase  rotation 
proportional  to  dj^.  Compensating  for  this  phase  rotation,  and  then 
averaging  the  short  term  spectral  estimates  yields  the  following  estimate 
of  the  cross  power  spectrum, 

Gr  r  ^  V  ^r  r 

12  fcf  12 

=  a  G^2(f)  .  (5-14b) 


84 


Multiplying  (5-14)  by  the  appropriate  weighting  function  and  computing  the 
IFT  of  the  product  results  in  an  estimate  of  the  GCC  function.  In 
practice,  d|^,  or  equivalently,  d'  =  b  -  1  is  not  known  and  it  is 
necessary  to  compute  (5-14)  for  several  hypothesized  values  of  d',  say 


X ' 1  ,  1  =  1,  2,  ...,  N,  then 


(5-15) 


r  (f) 


r  2 


Computing  the  IFT  of  (5-15)  for  1  =1,  ...,  N  results  in  an  ambiguity 
surface  which  is  a  function  of  the  hypothesized  values  for  d'  and  d^  (or 
B  and  D).  The  position  of  the  peak  on  this  surface  corresponds  to  the 
estimates  for  d'  and  d^.  This  technique  is  illustrated  in  Figure  5.2. 

5.4  Simulation  Results 

This  compensation  technique  was  implemented  and  preliminary 
simulation  results  comparing  the  SCC  and  SCOT  methods  have  been  obtained. 
Simulation  details  and  a  program  listing  is  given  in  Appendix  G.  Briefly, 
the  model  in  (5-2)  and  (5-3)  was  simulated  for  an  assumed  sampling 
frequency  of  2048  Hz,  T  =  4  seconds  =  8192  samples,  and  M'  =  33  segments 
(50%  overlap  was  used).  Each  segment  was  then  1/4  second  or  512  samples 
long.  The  initial  delay,  d^,  was  set  to  zero  at  the  start  of  the 
simulation  and  the  delay  rate,  d',  was  set  at  2.5  x  10“^  samples/sampling 
interval.  When  the  total  delay  became  greater  than  15  samples,  the  delay 
rate  was  set  to  -2.5  x  10“^  samples/sampling  interval  until  the  total 
delay  decreased  below  -15  samples  and  the  sign  of  d'  was  again  changed  and 
so  on.  Hypothesized  delay  rates,  x'-j,  of  0,  +2.5  x  10”^,  and  +5  x  10”^ 
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Figure  5.2  Frequency  Domain  Compensation  Technique 
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samples/sampling  interval  were  used.  The  SNR  was  -3  dB.  The  simulation 
results  for  the  SCOT  and  SCC  techniques  are  displayed  in  Figures  5.3  and 
5.4,  respectively.  In  both  cases.  Figure  (a)  shows  the  estimated  delay 
values  (circles  connected  by  dashed  line)  and  the  true  delay  values  (solid 
line)  as  a  function  of  time,  and  similarly.  Figure  (b)  shows  the  estimated 
delay  rate  (circles)  and  the  true  delay  rate  (solid  line)  as  a  function  of 
time.  Both  methods  provided  good  estimates  of  the  time  delay  with  the  SCOT 
predicting  the  delay  value  to  within  a  fraction  of  a  sample  in  every 
trial.  The  SCOT  method  did  a  much  better  job  of  estimating  the  delay  rate 
than  the  SCC  method.  This  is  probably  due  to  the  narrower  peak  of  the  SCOT 
function  compared  to  the  SCC  function,  however,  this  effect  has  not  been 
investigated  thoroughly.  Also  note  that  in  this  preliminary  simulation,  no 
interpolation  was  performed  to  obtain  the  estimates  of  the  delay  and  delay 
rate.  The  values  of  delay  and  delay  rate  corresponding  to  the  global  peak 
of  the  five  discrete  GCC  functions  obtained  from  the  processor  were  used  as 
the  estimates. 
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Figure  5.3  SCOT  Simulation  Results  -  Time-Varying  Delay 
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Figure  5.4  SCC  Simulation  Results  -  Time-Varying  Delay 
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CHAPTER  6 


SUMMARY  AND  RECOMMENDATIONS 

This  work  has  presented  a  detailed  study  of  time  delay  estimator 
performance  to  obtain  a  better  understanding  of  the  fundamental  limitations 
in  the  time  delay  estimation  problem.  To  this  end,  an  investigation  of  the 
behavior  of  the  Cramer-Rao  lower  bound  (CRLB)  and  the  recently  proposed 
correlator  performance  estimate  (CPE)  and  Ziv-Zakai  lower  bound  (ZZLB)  was 
conducted.  Derivations  of  these  three  performance  estimates  were  outlined 
and  expressions  for  the  variance  of  the  time  delay  estimate  for  specific 
signal  and  noise  power  spectra  were  obtained. 

It  was  observed  that,  although  the  ML  estimator  is  asymptotically 
efficient  in  attaining  CRLB  performance,  actual  performance  can  be  much 
worse  than  that  predicted  by  the  CRLB  for  a  given  SNR  and  observation 
time.  The  CPE  and  the  ZZLB  were  seen  to  be  characterized  by  a  threshold 
SNR.  Above  the  threshold  SNR,  the  CPE  and  the  ZZLB  coincide  with  the 
CRLB.  Below  the  threshold,  large  estimation  errors  or  anomalous  estimates 
become  dominant  and  the  performance  predicted  by  the  CPE  and  the  ZZLB 
degrades  rapidly  relative  to  the  CRLB.  As  SNR  goes  to  zero,  the  CPE  and 
the  ZZLB  are  further  characterized  by  a  prior  information  limit.  That  is, 
the  predicted  variance  approaches  a  constant  value  which  is  a  function  of 
the  maximum  possible  delay  value.  The  threshold  SNR  was  shown  to  be 
approximately  inversely  proportional  to  the  square  root  of  the  observation 
(coherent  processing)  time.  Simulation  results  were  found  to  be  in  good 
agreement  with  the  theoretical  predictions. 
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This  threshold  effect  was  found  to  have  important  implications 
related  to  coherent  versus  incoherent  processing  considerations.  A  simple 
incoherent  processor  was  considered  which  divides  the  T  second  observation 
interval  into  N  sections,  each  of  length  T/N  seconds.  Each  T/N  second 
section  is  coherently  processed  separately  to  obtain  N  time  delay  estimates 
and  finally  the  N  estimates  are  averaged  to  yield  a  single  estimate  at  the 
end  of  the  observation  interval.  Compared  to  a  coherent  processor,  which 
processes  the  entire  T  second  observation  interval  coherently,  the 
incoherent  processor  suffers  a  loss  of  5  log  N  dB  in  the  threshold  SNR, 
that  is,  (dB)  =  ~  ^  example,  if  the 

incoherent  processor  divides  100  seconds  of  data  into  100,  1  second 
sections,  the  coherent  processor  maintains  CRLB  performance  for  -  10  dB 
below  the  threshold  SNR  for  the  incoherent  processor.  Thus,  significant 
performance  gains  can  be  obtained  by  increasing  the  coherent  processing 
time. 


However,  increasing  the  coherent  processing  time  requires  a 
significantly  more  complex  structure  for  time  delay  estimators,  in  the  case 
of  a  time-varying  time  delay.  The  ML  estimator  for  the  time-varying  delay 
case  was  seen  to  require  that  the  received  signals  be  pre-processed  to 
compensate  for  the  relative  time  compression  between  the  two  signals.  In 
general,  the  required  pre-processing  is  difficult  to  achieve  from  an 
implementation  viewpoint,  in  addition  to  being  computationally  intensive. 
Therefore,  a  simplified  compensation  scheme  was  proposed  which  essentially 
computes  short-term  correlation  functions  so  that  the  time  delay  can  be 
considered  constant  over  the  correlation  interval.  The  short-term 
correlation  functions  are  then  combined  coherently  to  achieve  a  longer 
coherent  processing  time.  Time  domain  (actually,  time  delay  domain) 
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compensation  techniques  which  coherently  sum  the  short-term  correlation 
functions  directly  are  discussed  in  [49,50].  The  method  described  in  this 
work  is  implemented  in  the  frequency  domain  which  is  attractive  in  that  it 
is  easily  incorporated  into  the  structure  of  the  GCC  processor. 

Preliminary  simulation  results  demonstrate  that  this  compensation  method  is 
effective  and  allows  the  delay  and  delay  rate  to  be  accurately  estimated. 
However,  the  method  has  only  been  tested  for  a  limited  set  of  conditions 
and  further  testing  is  required. 

Future  work  based  on  this  dissertation  could  take  many  directions. 
Some  of  these  include: 

1)  Extension  of  the  time  delay  model  in  (2-1)  to  account  for 
multiple  sources  or  multi-path  interference  and  the  development 
of  the  corresponding  expression  for  the  ZZLB.  (The  CRLB  for  the 
multi-sensor,  multi-target  problem  is  developed  by  Ng  in  [51].) 

2)  More  extensive  testing  of  the  proposed  compensation  technique  for 
the  time-varying  time  delay  model  and  comparison  with  time  domain 
implementations. 

3)  Consideration  of  the  relative  importance  of  the  fundamental 
limitations  on  time  delay  estimator  performance  discussed  here 
(e.g.,  due  to  SNR  and  observation  time)  with  the  limitations 
imposed  by  uncertainty  in  sensor  position. 

These  and  other  problems  of  equal  importance  provide  the  basis  for  a  great 
deal  of  exciting  research  yet  to  be  done  in  time  delay  estimation. 
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APPENDIX  A 


THE  ML  ESTIMATOR-FIXED  TIME  DELAY 

The  ML  estimator  for  the  case  of  a  stationary  or  fixed  time  delay  is 
derived  by  Knapp  and  Carter  in  [6].  The  derivation  is  summarized  in  this 
appendix  for  convenience.  The  mathematical  model  is  that  of  Section  2.1 
and  the  same  assumptions  are  made. 

Suppose  an  observation  vector,  which  is  dependent  upon  the 
parameter  to  be  estimated,  say  x,  is  given.  For  the  TDE  problem,  the 
parameter  of  interest  is  the  true  time  delay  D.  The  ML  estimate,  is 

the  value  of  x  which  is  most  likely  to  have  resulted  in  the  occurrence  of 
the  observed  vector,  More  precisely,  the  ML  estimate  is  given  by 
Dml  =  Xq  such  that 

I  ■‘^0^  —  I  allowable  x,  (A-1) 

where  p(^  |  x)  is  the  conditional  probability  density  function  of  ^  given  x. 

For  purposes  of  analysis,  consider  periodic  extensions  of  the 
received  signals,  r^(t)  and  r2(t),  with  period  T,  where  T  also 
represents  the  observation  time.  Then,  for  any  integer  m 

r.(t  +  mT)  =  r^-(t)  ,  0  <  t  <  T,  i  =  1,2.  (A-2) 

The  Fourier  series  coefficients  of  the  received  signals  can  be  computed  as 
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( A-3a) 


T 


~  1  S  '"i 
0 


where 


=  2  IT  /  T . 


(A-3b) 


The  time  waveforms  can  be  reconstructed  from  the  Fourier  coefficients  using 
the  relationship 


CO 


(A-4a) 


=  r . (t)  ,  0  <  t  <  T  . 


(A-4b) 


In  practice,  the  received  signals  are  band  limited  so  the  number  of 
non-zero  Fourier  coefficients  will  be  finite,  say  -N  £  k  <  N.  Thus,  as 
seen  in  (A-4)„  r^.  (t)  is  characterized  by  the  set  of  Fourier  series 
coefficients  <R^(k),  -N  £  k  £  nI  . 


Defining  the  vector  R(k)  =  [R2^(k),  R2(k)]',  where  '  denotes 
transpose,  the  ML  estimate  of  the  time  delay  is  obtained  by  maximizing  the 
joint  density  function,  p(^(-N),  ^(-N+1),  ...,  ^(N)  |  x),  with  respect  to 
the  time  delay  variable  x.  Observe  that  R^-(k)  is  a  linear  transformation 
of  '"^■(t),  which  is  Gaussian,  so  R^-(k)  is  also  Gaussian  for  i  =  1,2  and 
-N  <  k  <  N.  Further,  it  can  be  shown  that  [36,  p.  461] 


E[Rl(k)  R^d)]  = 


(A-5) 


0 


k  ^  1  . 
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It  is  also  easily  seen  that  E[R^.(k)]  =  0  ,  since  r^-(t)  is  zero  mean. 
Therefore,  the  vectors,  ^(k),  -N  £  k  £  N,  are  uncorrelated,  zero  mean, 
Gaussian  random  vectors.  Let  ^  =  (£(-N),  ...,  £(N)),  then  the  joint 
density  function  can  be  written  as 


N 


p(R  I  t)  =  I  >  p(R(k)  I  t) 
k=-N 


(A-6a) 


where 


p(R(k)  I  t)  =  [2t[ 


exp 


-  7  a*'(k)£j;i\a(i<)j .  (A-6b) 


In  (A-6b)  ,^1^1  ^  denotes  the  covariance  matrix  of  the  vector  £(k)  given  t, 
and  is  defined  as 


EkI,  =  E  [R(k)  R*'(k)]T] 


G*  ( kw  ) 

(kWp), 

(A-7a) 


(A-7b) 


(A-7c) 


where  Q  is  defined  to  be  a  spectral  density  matrix.  The  power  spectra  in 
(A-7b)  are  assumed  known.  The  cross  power  spectra  are  dependent  on  the 
hypothesized  time  delay  t,  hence, is  a  function  of  t.  The 
determinant  of^j^j^,  dropping  the  frequency  argument  of  the  power 
spectra,  is 
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Skit 


=  (G  +  G  )(a^G  +  G  ) 

'  ss  '^1^1  ^2^2 


-jkw  T  jkw  T 

-  (“^55  6  ^  ^“^SS  ®  ^  ’ 


(A-8) 


which  is  seen  to  be  independent  of  the  hypothesized  delay  Substituting 
(A-6b)  into  (A-6a)  yields 


% 


P(R 


t)  =  c  exp  ( 


-1 

T 


Jl) 


where 


and 


k=-N 


R*’(k)E,^ 


-1 


,  i(k) 


N 


c  = 

k=-N 


(A-8a) 


(A-8b) 


(A-8c) 


Since  c  does  not  depend  on  t,  maximizing  p(^|t)  with  respect  to  t  is 
equivalent  to  minimizing  since  E[Jj^]  >  0. 


At  this  point,  an  approximation  is  made  based  on  the  relationship 
between  the  Fourier  series  and  the  Fourier  transform  (FT).  For  large  T,  it 
can  be  shown  that  [52,  pp.  23-25]. 


cx> 

=  y  R*'(f)  Q"^(f)  R(f)  df 

CO 


(A-9) 
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where  R(f)  =  [R]^(f),  R2(f)]'  and  where  R^(f)  represents  the  FT  of 
r^(t).  Evaluating  Q~^(f)  from  (A-7)  and  substituting  into  (A-9), 
can  be  written  as 


Jl  -  J2  -  J3 


(A-lOa) 


where 


^2- 


Rl(f) 


^1^1 


TT 


R2(1^) 


^2^2 


2  1 

TT 


r^^TTTTTT 

rir2 


df  , 


(A- 10b) 


/*  *  «  2  (f)  e- 

'^3  =  V  ^  (f)  ^  (f)  (1  -  C,  ,  (f)) 


j2TrfT 


’"2'"2 


'"l'"2 


and  recall 


^r  r  =  G 
rir2  b 


Gr  r  (f) 

'^r  2 


-  -  (TTGrTTf) 

'^r  1  '^2'^2 


(aGss(f)) 

G.  .  (f)  G  (f) 

'^r  1  '^2'^2 


(A-lOd) 


Now  does  not  depend  on  t,  therefore,  maximizing  J3  with  respect  to  t 
is  equivalent  to  minimizing  Jp  and  yields  the  ML  estimate.  The 
information  about  the  true  delay,  D,  is  contained  in  R]^(f)  and  R2(f) 
which  are  obtained  from  the  received  signals,  r^(t)  and  r2(t),  over  the 
observation  interval  T.  The  other  terms  in  (A-lOc)  depend  on  the  known 
power  spectra  with  the  hypothesized  delay  t  appearing  in  the  exponential 
term.  The  product  R2^(f)  R2(f)  can  be  viewed  as  T  times  an  estimate 
of  the  cross  power  spectrum,  G  (f). 

'^r  2 

Re-writing  in  terms  of  G  (f)  and  C  (f)  and  noting  that 
j  2  ^r?. 

multiplying  by  a  positive  constant  does  not  affect  the  location  of  the 
peak,  the  ML  estimate  of  the  time  delay  is  the  value  of  t  at  which 
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G.  „  (f) 
2 


2 


'^l'^2 


m 


1  -•  c 


if) 

m 


(A-11) 


attains  its  peak  value.  From  (A-11),  it  is  seen  that  the  ML  weighting 
function  is  indeed  given  by  (2-19), 


C.  .  (f) 


W(f)  = 


1  ''l''2 


G.  .  (f)  •  1  -  (f) 


’'1^2 


(A-12) 
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APPENDIX  B 


SIMULATION  DETAILS-FIXED  TIME  DELAY 


A  computer  simulation  was  conducted  to  corroborate  the  theoretical 
TDE  performance  predictions  and  to  compare  the  performance  of  several  of 
the  GCC  processors.  The  simulation  results  are  presented  in  Section  4.3. 
The  simulation  was  implemented  on  a  VAX  11-750. 

The  signal  model  of  (2-1)  was  simulated  as  follows.  Two  independent, 
Gaussian,  pseudo-random  sequences  having  zero  mean  and  unit  variance  were 
generated  to  represent  the  additive  noise  sequences,  n^(t)  and  n2(t). 

A  third  independent  Gaussian  sequence,  again  with  zero  mean  and  unit 
variance,  was  generated  to  represent  the  broadband  source  signal,  s(t). 

The  required  Gaussian  random  sequences  were  obtained  by  first  generating 
random  sequences  which  were  uniformly  distributed  on  the  interval  (0,  1) 
using  the  technique  in  [53,  p.  S-11].  These  uniformly  distributed  random 
sequences  were  then  transformed  to  Gaussian  random  variables  using  the 
method  in  [54,  p.  953].  Specifically,  if  and  U2  represent  a  pair  of 
random  variables  uniformly  distributed  on  (0,  1),  the  transformation 


X^  =  (-2  In  cos(2itU2) 


(B-la) 


X2  =  (-2  In  sin(2iTU2) 


(B-lb) 


yield  a  pair  (Xp  X2)  of  independent  Gaussian  random  variables  with 
zero  mean  and  unit  variance. 
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The  three  Gaussian  random  sequences  were  then  processed  through  a 
12th  order  (six  cascaded  2nd  order  sections)  low-pass  Butterworth  filter 
with  a  cut-off  frequency  of  100  Hz,  relative  to  an  assumed  sampling 
frequency  of  2048  Hz.  The  noise  sequences  were  scaled  to  obtain  the 
desired  SNR.  To  obtain  r^(t),  the  signal  sequence  was  simply  added  to 
one  of  the  noise  sequences,  and  to  obtain  r^it),  the  signal  sequence  was 
delayed  (0=4  samples)  and  added  to  the  second  noise  sequence. 

The  GCC  function  of  the  sequences,  r^(t)  and  r2(t)  was  computed 
following  the  procedure  given  in  [53,  pp.  2.3-1  -  2.3-18].  The  sequences 
were  processed  using  512  point  data  segments  and  50%  overlap.  Thus,  each 
data  segment  represents  1/4  second  of  data.  To  obtain  a  2  second 
observation  time,  17  overlapped  segments  must  be  processed.  Similarly,  an 
8  second  observation  time  requires  65  overlapped  segments.  Each  512  point 
segment  was  weighted  with  a  Hamming  window  and  then  zero  padded  to  1024 
points.  The  Fourier  transform  of  each  1024  point  segment,  for  both  r^(t) 
and  r2(t),  was  computed  using  a  fast  Fourier  transform  (FFT).  The 
resulting  Fourier  coefficients  were  used  to  compute  auto  and  cross-power 
spectral  estimates  for  each  segment,  and  these  estimates  were  then  averaged 
to  obtain  the  final  estimates  of  the  spectra  for  the  observation  interval. 

Using  the  spectral  estimates,  the  GCC  weighting  functions  were 
computed  and  multiplied  by  the  estimate  of  the  cross  power  spectrum.  The 
inverse  FFT  of  the  result  was  computed  to  yield  an  estimate  of  the  GCC 
functions.  These  functions  were  searched  for  their  peak  values  for  delay 
values  over  the  range  +128  samples.  Thus,  the  value  of  in  seconds  is 
128/2048  =  1/16  second.  The  time  delay  estimate  is  given  by  the  delay 
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value  corresponding  to  the  peak  of  the  GCC  function  in  the  search 
interval.  The  simulation  was  conducted  for  two  observation  times,  2  and  8 
seconds,  and  over  a  range  of  SNR  values  from  0  to  -15  dB.  A  total  of  2000 
trials  was  conducted  at  each  SNR  to  obtain  the  experimental  time  delay 
variances,  which  are  displayed  in  Figures  4. 4-4. 8. 
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APPENDIX  C 


NUMERICAL  EVALUATION  OF  P[A] 

In  this  appendix,  techniques  for  numerically  evaluating  the 
expression  for  P[A]  in  (3-16a)  are  discussed.  Recall 


PEA]  =  1 


e-(x-3r)^/2 


(C-1) 


where  x  =  f(B,  T,  SNR)  and  x  =  g(SNR)  are  given  in  (3-16).  This  equation 
can  be  re-written  in  terms  of  the  Q  function  to  give 


P[A]  =  1  - 


/ 

*CD 


1  g-(x-x)^/2 

N/2tt 


CQ(-xx)]'^  ^  dx 


(C-2a) 


where 

CD 

Q(x)  =  — ^  / e~^  dy  , 


(C-2b) 


and  note  that  Q(-<x>)  =  1,  Q(0)  =  UZ,  and  Q(-x)  =  1  -  Q(x). 

Consider  the  integral  on  the  right  hand  side  of  (C-2a).  Making  the 

change  of  variable  u  =  (x-x)/  2  yields 

oo  2 

1  -  P[A]  =  —  f  Q(-x(2u  +  x))^“^  du.  (C-3) 

L 

This  expression  for  1  -  P[A]  can  be  approximated  using  the 
Gauss-Hermite  Quadrature  formula  [54,  p.  890,  eq.  25.4.46].  Consider  a 
function  6(x) ,  then 
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e~^  G(x)  dx  = 

That  is,  the  integral  can  be  approximated  by  a  weighted  sum  of  the  function 

G(x)  evaluated  for  specific  arguments.  In  (C-4),  x^  ^  represents  the  ith 

zero  of  the  mth  order  Hermite  polynomial,  and  W  -is  the  weight 

m,  I  ^ 

corresponding  to  x^  .  and  is  a  function  of  the  mth  order  Hermite 
polynomial.  If  the  zeroes  are  arranged  in  ascending  order  and  the  weights 
are  arranged  in  corresponding  order,  then  the  following  relations  hold. 


m,  1 


G(x  . ) 
m,  1 


(C-4) 


/ 


~  ~  ^m,m+l-i  ’ 


and 


'^m,i  ^'^m, 111+1-1 


(C-5a) 

(C-5b) 


The  zeroes  and  weights  are  available  in  tables,  for  example  see  [55,  Table 
5].  Note  that  most  tables  make  use  of  the  property  of  (C-5)  and  only  list 
the  positive  zeroes  and  their  corresponding  weights. 


Applying  the  approximation  of  (C-4)  to  (C-3)  yields 

1  .  P[A]  =  i  [Q(-x(^  x^_,  -  .  (C-6) 

At  the  end  of  this  appendix,  a  listing  is  given  of  the  FORTRAN 
program  BOUND. FOR,  which  computes  the  CRLB,  the  CPE,  and  the  ZZLB  as 
functions  of  SNR.  Additionally,  listings  of  the  subroutines  PROB.FOR  and  Q. 
FOR,  which  evaluate  P[A]  and  the  Q  function,  respectively,  are  included. 

The  values  of  P[A]  which  are  of  greatest  interest  are  of  the  order  10“^ 
or  smaller.  Thus,  the  approximation  in  (C-6)  must  be  evaluated  to  5  or 
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more  significant  digits.  This  requires  a  large  value  for  m  (m  =  60  in  the 
subroutine)  and  the  use  of  double  precision  calculations.  An  alternative 
method  of  evaluating  P[A]  is  proposed  by  Nuttall  in  [45],  which  reduces  the 
number  of  significant  digits  required  for  the  calculation.  These  two 
techniques  have  been  found  to  yield  very  similar  values  for  P[A]. 
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C 

C  CRLB  -  CPE  -  ZZLB  ROUTINE 

C 
C 
C 
C 
C 
C 
C 
C 
C 


VAX-11  FORTRAN  SOURCE  FILENAME: 


DEPARTMENT  OF  ELECTRICAL  ENGINEERING 


BOUND. FOR 


KANSAS  STATE  UNIVERSITY 


REVISION 


DATE 


PROGRAMMER(S) 


00.0 
01  .0 


SCARBROUGH 

SCARBROUGH 


JUNE  1983 
JULY  1983 

C 


PURPOSE 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

Q*********************************************************************** 

C 


This  program  computes  the  Cramer-Rao  lower  bound  (CRLB), 
lanniello's  estimate  of  Correlator  Performance  (CPE), 
and  the  Ziv-Zakai  lower  bound  (ZZLB),  for  the  case  of 
low  passed  gauss  Ian  signal  &  noise  sequences. 

ROUTINE (S)  ACCESSED  OR  CALLED  BY  TH I S  ROUTINE 
PROB,  0,  header,  SGOPEN,  SGTRAN 


INPUT  PARAMETERS 
SF 

FREQ  = 
TIME  = 
WIDTH  = 
SNR1  = 
SNR2  = 
SSTEP  = 
EPS 

PSTEP  = 


SAMPLING  FREQUENCY 
CUTOFF  FREQUENCY 
OBSERVATION  TIME 
CORRELATION  WINDOW  WIDTH 
MIN  SNR  VALUE (DB) 

MAX  SNR  VALUE (DB) 

SNR  STEP  SIZE 
EPSILON  VALUE 

STEP  SIZE  FOR  PRINTING  CRB/ CPE  VALUES 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

Q******************************1F**************************************** 

C 

C  NOTE  :  The  equations  for  the  CPE  and  the  ZZLB  assume 

C  that  the  signal  has  a  low-passed  (flat)  power  spectrum. 

C 

c 

REAL  CRB(8002),  CPE(8002),  ZZB(8002),  SNRSAV(4001 ) ,  PASAV(4001) 

REAL  BW,  CPELIM,  CRBVAL,  ZZBLIM,  FREQ,  EPS 

REAL  ETA2,  ETAZ,  ETAZ2,  PA,  PE,  Q_ETAZ 

REAL  SNR,  SNR1,  SNR2,  SSTEP,  SNRDB,  SNRFCT,  SQSNR1 

REAL  TIME,  ZBW,  ZBWRD,  ZZ1 ,  ZZ2 

INTEGER  MEXP,  NPT,  NPT2,  PSTEP,  WIDTH,  UN  I  TP 

LOGICAL  EQUAL1,  EQUAL2 

CHARACTER* 40  CRBNAM,  CPENAM,  ZZBNAM 

DATA  PI  /3. 14592654/ 
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c 

C  INITIALIZE  CONSTANTS 
C 

NPT  =  0 

EQUAL  1  =  .FALSE. 

EQUAL2  =  .FALSE. 

TWOPI  =  2.  *  PI 
SQTPI  =  SQRT( TWOPI) 

UNITP  =  99 
C 

C*******************************************^**^*^**^***^***^* 


c 

C  FORMAT  STATEMENTS 
C 


1  FORMAT(A) 

3  FORMATS  ’,A,F10.2,A) 

5  FORMATS  ’,A,I10,A) 

7  FORMATS  ’,A,A) 

9  FORMATC  » ,/ ,4X, » SNR( dB) » ,9X, »CRLB » , 1 2X, »CPE M 1 X, »ZZLB ’ ,/) 

11  FORMAT( *1 ',//,4X, ’SNR(dB) ’,6X, *PR(anomaly) ’,4X, *ZZB  Term  1', 

&  5X,»ZZB  Term  2’,/) 


(2*******************************************ih^**********^***** 


C 

C 

c 


c 


ACCEPT  INPUT  PARAMETERS 


TYPE  1,  *$SAMPLING  FREQUENCY 
READ  *,  FS 

TYPE  1,  ’SOBSERVATION  TIME (sec) 
READ  *,  TIME 

TYPE  1,  'SCORR  W I NDOW( samples) 
READ  *,  WIDTH 
TYPE  1,  ’SCUTOFF  FREQ 
READ  FREQ 

TYPE  1,  ’$MIN  SNR  VALUE (dB) 

READ  *,  SNR1 

TYPE  1,  ’$MAX  SNR  VALUE (dB) 

READ  *,  SNR2 

TYPE  1,  ’$SNR  STEPSIZE(dB) 

READ  SSTEP 

TYPE  1,  ’SEPSILON  VALUE 

READ  *,  EPS 

TYPE  1,  'SPRINT  STEPS  I ZE 
READ  *,  PSTEP 


C 

C  COMPUTE  PROGRAM  VARIBLES 
C 


BW  =  FREQ 
ZBW  =  BW  *  2. 

DO  =  WIDTH  /  FS 
CPELIM  =  (D0**2)  /  12. 

ZZBLIM  =  2.  *  CPELIM 
MEXP  =  IFIX(2.  *  DO  *  BW) 

CRBVAL  =  3.0  /  (8.0  *  (Pl*«2)  *  (BW**3)  *TIME) 
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ooo  oooo  o  o  o  ooo  oooooo  o  ooo 


COMPUTE  CRLB,  CPE  &  ZZLB  AS  FUNCTIONS  OF  SNR 

DO  SNRDB  =  SNR1,  SNR2,  SSTEP 
NPT  =  NPT  +  1 
JPT  =  2  *  NPT  -  1 
IPT  =  2  *  NPT 

SNRSAV(NPT)  =  SNRDB 
SNR  =  10.**(SNRDB  /  10) 

SNRFCT  =  (SNR**2)  /  (1.  +2.  *  SNR) 

SQSNR1  =  SQRTCSNRFCT  +1.) 

COMPUTE  CRLB 

CRBTMP  =  CRBVAL  /  SNRFCT 
COMPUTE  CPE 

PA  =  PROB(SNR,BW,TIME,MEXP) 

PASAV(NPT)  =  PA 

CPETMP  =  PA  *  CPELIM  +  (1.  -  PA)  *  CRBTMP 
COMPUTE  ZZLB 

ZTMP1  =  LOG((1.  +  SQSNRl)  /  2.) 

ZTMP2  =  (SQSNRl  -  1.)  /  (2.  *  SQSNRl) 

ZTMP3  =  EXP(-ZBW  *  TIME  *  (ZTMP1  -  ZTMP2)) 

ETAZ  =  SQRTCZBW  *  TIME  *  (SQSNRl  -  1.)  /  SQSNRl) 
ETAZ2  =  ETAZ**2 
Q_ETAZ  =  Q(ETAZ) 

ZTMPA  =  (ETAZ2  -  1)  *  Q_ETAZ 

ZTMPB  =  ETAZ  *  EXP(-ETAZ2  /  2.)  /  SQTPI 

ZTMPC  =0.50 

PE  =  Q_ETAZ  *  ZTMP3 
ZZl  =  PE  *  ZZBL IM  '  ■ 

ZZ2  =  2  *  CRBTMP  *  (ZTMPA  -  ZTMPB  +  ZTMPC) 

ZZBTMP  =  ZZl  +  ZZ2 

COMPUTE  LOGIO  OF  STANDARD  DEVIATION  OF  TIME  DELAY  ESTIMATE 

CRBTMP  =  LOGIO (CRBTMP)  /  2. 

CPETMP  =  LOGIO (CPETMP)  /  2. 

ZZBTMP  =  LOGIO (ZZBTMP)  /  2. 

DETERMINE  TWRESHOLD  SNR  VALUES 

IF  (. NOT. EQUAL! )  THEN 

DIFF  =  CPETMP  -  CRBTMP 
IF  (DIFF  .LT.  ABS(EPS  *  CRBTMP))  TMEN 
EQUAL!  =  .TRUE. 

INDEX!  =  NPT 
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END 


END 

IF 


INDEX2  = 
IF 


IPT 


IF  ( 


.N0T.EQUAL2)  THEN 
DIFF  =  ZZBTMP  -  CRBTMP 


END 


IF 


END 

IF 


(DIFF  .LT. 
EQUAL2  = 
INDEX3  = 
INDEX4  = 
IF 


ABS(EPS 

.TRUE. 

NPT 

IPT 


*  CRBT>1P))  THEN 


C  STORE 
C  (SNR 
C 


SNR  VALUES  &  BOUND  VALUES  IN  OUTPUT  ARRAY 
VALUES  NEEDED  FOR  XY  PLOTTING) 


CRB( JPT) 
CRB( IPT) 
CPE (JPT) 
CPE( IPT) 
ZZB( JPT) 
ZZB( IPT) 


SNRSAV(NPT) 

CRBTMP 

SNRSAV(NPT) 

CPETMP 

SNRSAV(NPT) 

ZZBTMP 


END  DO 


',CRBNAM,'REAL',NPT2) 
’,CPENAM, 'REAL’,NPT2) 
»,ZZBNAM, 'REAL’,NPT2) 


C 

C  OPEN  OUTPUT  FILES,  OUTPUT  CRLB  &  CPE  VALUES  TO  DISK 
C 

NPT2  =  NPT  *  2 

CALL  SGOPEN(l,»WRITE’,’CRB  FILE 
CALL  SG0PEN(2, »WRITE»,’CPE  FILE 
CALL  SG0PEN(3, 'WRITE*, 'ZZB  FILE 
CA LL  SGTRAN ( 1 , ' WR I TE ' , ' RE AL ' , CRB , NPT2 ) 

CALL  SGTRAN(2, 'WRITE' , 'REAL' , CPE, NPT2) 

CALL  SGTRAN(3, 'WRITE' , 'REAL' ,ZZB,NPT2) 

C 

C 

C  CREATE  LISTING  FILE:  WRITE  HEADER 
C  ECHO  INPUT  PARAMETERS 

C 
C 
C 
C 


OUTPUT  CRB  -  CPE  INTERCEPT  VALUES 
OUTPUT  CRB  -  ZZB  INTERCEPT  VALUES 
OUTPUT  CRB,  CPE  &  ZZB  VALUES 

CALL  HEADER  (UNITP,'***  CRLB  -  CPE  -  ZZLB  PROGRAM 

***  1 ) 

WRITE(UNiTP,3) 

*  Sampling  Frequency 

f 

P 

FS  , 

'  Hz' 

WRITE(UNITP,3) 

'  Observation  Time 

» 

P 

TIME  , 

'  sec* 

WRITE(UNITP,5) 

*  Correlation  Window 

» 

P 

WIDTH, 

'  samp 

WRITE(UNITP,5) 

'  Ind  Corr  Values 

f 

P 

MEXP 

WRITE(UNITP,3) 

'  Cutoff  Freq 

» 

P 

FREQ  , 

'  Hz* 

WRITE(UNITP,3) 

'  Min  SNR  Value 

f 

P 

SNR1  , 

*  dB* 

WRITE(UNITP,3) 

'  Max  SNR  Value 

f 

P 

SNR2  , 

'  dB* 

WRITE(UNITP,3) 

'  SNR  Steps ize 

f 

P 

SSTEP, 

*  dB* 

WRITE(UNITP,*) 

*  Eps i 1  on  Value 

f 

P 

EPS 

I  es ' 


WRITE(UNITP,*) 
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WRITE(UNITP,3)  ' 

Threshold  SNR(dB) 

.  1 
•  P 

SNRSAVC INDEX1 ), 

WRITE(UNITP,*)  ' 

Prob  of  Anomaly 

.  1 
•  P 

PASAVC INDEX 1 ) 

WRITE(UNITP,*)  ' 

CPE  Value 

.  1 
•  P 

CPE( INDEX2) 

WRITECUNITP,*)  ' 

CRLB  Value 

.  1 
•  P 

CRB( INDEX2) 

WRITE(UNITP,*) 
WRITE(UNITP,3)  ' 

Threshold  SNR(dB) 

.  1 
•  P 

SNRSAVC INDEX3), 

WRITE(UNITP,*)  ' 

CPE  Value 

.  1 
•  P 

CPEC INDEX4) 

WRITE(UNITP,*)  ' 

ZZLB  Value 

•  1 
•  P 

CRB( INDEX4) 

WRITE(UNITP,*) 
WRITE(UNITP,7)  ' 

CRLB  Filename 

.  1 
•  P 

CRBNAM 

WRITE(UNITP,7)  ' 

CPE  Filename 

.  1 
•  P 

CPENAM 

WRITE(UNITP,7)  ' 

ZZLB  Fi 1 ename 

.  1 
•  P 

ZZBNAM 

WRITE(UNITP,9) 

DO  IPT  =  1,  NPT2, 

PSTEP*2 

JPT  =  IPT  +  1 

WRITE(UNITP,*) 

CRB(IPT),  CRB(JPT), 

CPE (JPT),  ZZB(JPT) 

END  DO 
C 

STOP 

END 


'  dB' 


'  dB* 
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c*** ********************************************************* *********** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


PROBABILITY  OF  ANOMALY  ROUTINE 
VAX-11  FORTRAN  SOURCE  FILENAME: 
DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
REVISION  DATE 

00.0  MAY  1983 


PROB.FOR 

KANSAS  STATE  UNIVERSITY 
PROGRAMMER(S) 

KENT  SCARBROUGH 


********************************************************************** 


C 

C  CALLING  SEQUENCE 

C 

C  VALUE  =  PROB(SNR,BW,T,M) 

C 

C  PURPOSE 

C 

C  This  routine  calculates  the  probability  of  anomaly 

C  for  band-limited  signal  and  noise  sequences.  The 

C  signai  &  noise  are  assumed  to  have  the  same  ideai 

C  band( low)-passed  spectrai  characteristics. 

C 

C  ROUTINE (S)  ACCESSED  OR  CALLED  BY  THIS  ROUTINE 

C 

C  Q 

C 

C  ARGUMENT(S)  REQUIRED  FROM  THE  CALLING  ROUTINE 

C 


c 

p 

SNR 

SIgna l-to-Noi se  ratio 

L/ 

c 

n 

BW 

Bandwidth  of  signal  and  noise 

c 

p 

T 

Observation  time  (seconds) 

c 

M 

Number  of  independent  correlation  estimates 

C 

C  ARGUMENT(S)  SUPPLIED  TO'  THE  CALL ING  ROUTI NE 

C 

C  NONE 

C 

0*********************************** ****************************** ****** 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NOTE  1:  The  equation  for  PR(anomaly)  was  derived  following 
lannieilo's  paper  in  IEEE  TRANS  ASSP,  Dec  1982. 

NOTE  2:  The  formula  for  the  Gauss-Hermi te  Quadrature  rule 

can  be  found  In  "Handbook  of  Mathematical  Functions" 
by  Abramowitz  &  Stegun,  Nat' I  Bureau  of  Standards,  1964. 
See  equation  25.4.46,  page  890. 

NOTE  3:  The  weights  and  zeroes  for  the  G-H  Quad,  rule  were 
taken  from  "Gaussian  Quadrature  Formulas"  by  Stroud 
&  Secrest,  Prentice  Hall,  1966.  See  Table  5. 


no 


C 

REAL  FUNCTION  PROB(SNR,BW,T,M) 

C 

DOUBLE  PRECISION  WEIGHT(60),  ZERO(60) 

LOGICAL  FIRST/. TRUE./ 

INTEGER  M,  I,  IWT,  JWT,  NWT 
REAL  ALPHA,  BETA,  BW,  P,  PBAR,  PTMP 
REAL  SNR,  T,  X 
C 

DATA  SQRTP I /I .772453851/,  SQRT2/1 .41 421 3562/ 

C 

Q**** ******************************************************  ************* 

C 

DATA  (ZERO( I ), 1=1 ,30)/  10.15910924618,9.520903677013, 

&  8.992398001405,8.5205692841 18,8.085188654249, 

&  7.675839937505,7.286276594396,6.912381532189, 

&  6.551 259167063,6.200773557993,5.859290196394, 

&  5.525521086139,5.198426534576,4.877150077473, 

&  4 . 560973757936 , 4 . 249286435956,3. 941 560733926 , 

&  3.637335876171 ,3.336204653548,3.037803338231 , 

&  2.741803748070,2.447906902308,2.155837871230, 

&  1 .865341531233,1 .57617901 1 975, 1 .2881 24674869, 

&  1  .000963499561 ,0.714488781673,0.428500064221 , 

&  0.142801238703  / 

C 

Q*********************************************************************** 


C 

C  NOTE:  THE  FIRST  TWO  (&  LAST  TWO)  WEIGHTS  HAVE  BEEN  SET  TO 
C  ZERO  DUE  TO  THE  DYNAMIC  RANGE  LIMITATIONS  (E-JB  to  E+38) 

C  OF  FORTRAN  77  ON  THE  VAX  WITHOUT  THE  FLOATING_G  OPTION. 

C  THE  ACTUAL  VALUE  OF  THE  WEIGHTS  SHOULD  BE: 

C 

C  WEIGHT(I)  =  .1 109587247 97E-44  =  WEIGHT(60) 

C  WEIGHT(2)  =  . 24397475881 5E-39  =  WEIGHT(59) 

C 

C 

DATA  (WEIGHT( I ), 1=1 ,30)  /  .OOOOOOOOOOOOE-00, . OOOOOOOOOOOOE-00, 
&  .37716267271 2E-35 ,. 1 33255961 1 76E-31 , . 1 71 55731 4767E-28, 

&  .  1  0294059971 7E-25 , . 334575695575E-23 , . 651 256725750E-21 , 

&  .81 5364047302E-1 9, . 692324790958E-1 7 , .41 524441 0969E-1 5 , 

&  . 1 81 662457626E-1 3, . 594843051 606E-1 2, . 1 48895734906E-1 0 , 

&  .289935901 2808E-9 , . 4456822775226E-8, . 547555461 9277E-7, 

&  . 543351 61 34205E-6, . 4394286936267E-5 , . 291 8741 9041 56E-4, 

&  . 1 60277334681 8E-3 , . 731 7735569655E-3, . 2791 324828953E-2, 

&  . 89321 78360308E-2, . 24061 2727661  IE-1,. 5471 89709321 8E-1 , 

&  .1052987636977856,. 1717761569188851 ,.2378689049586589, 

&  .2798531175228290  / 

C 

Q*********************************************************************** 


C 

C 


SAVE  SQRT2,  SQRTP I,  FIRST,  WEIGHT,  ZERO 


NWT  =  60 
C 


C 
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ooooo  ooo 


C  ON  FIRST  CALL  TO  ROUTINE,  INITIALIZE  SECOND  HALF  OF 
C  WEIGHTING  AND  ZERO  ARRAYS 
C 

IF  (FIRST)  THEN 

DO  IWT  =  1,  NWT  /  2 

JWT  =  NWT  -  IWT  +  1 
WEIGHT(JWT)  =  WEIGHT(IWT) 

ZERO(JWT)  =  -  ZERO( IWT) 

END  DO 
END  IF 
C 


COMPUTE  CONSTANTS  -  FUNCTIONS  OF  SNR 

ALPHA  =  SQRT2  *  SQRTCBW  *  T)  *  SNR 
&  /  SQRT(  SNR*^f2  +  (1  +  SNR)^f^f2  ) 

BETA  =  SQRT(  1  +  (SNR**2)  /  (1  +  SNR)**2  ) 


USE  GAUSS-HERMITE  (QUADRATURE  TO  COMPUTE  P( anomaly) 


PBAR  =0.0 
DO  IWT  =  1,  NWT 

X  =  BETA  *  (SQRT2  ZERO(IWT)  +  ALPHA) 

P  =  1 .0  -  Q(X) 

PTMP  =  WEIGHT(IWT)  *  (P**(M-1))  /  SQRTPI 
PBAR  =  PBAR  +  PTMP 
•  END  DO 

PROB  =1.0-  PBAR 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Q  FUNCTION 

VAX- 11  FORTRAN  SOURCE  FILENAME: 
DEPARTMENT  OF  ELECTRICAL  ENGINEERING 


(UTILITY  LIBRARY) 

Q.FOR 

KANSAS  STATE  UNIVERSITY 


REVISION  DATE 


PROGRAMMER(S) 


00.0  MAR  31,  1983  F  W  RATCL I FFE 

01.0  APR  08,  1983  K  SCARBROUGH 


0*********************************************************************** 


C 

C  CALLING  SEQUENCE 

C 

C  VALUE  =  Q  (X) 

C 

C  PURPOSE 

C 

C  This  routine  evaluates  the  function  Q(x),  that  is, 

C  the  integral  fHom  x  to  +  infinity  of  a  Gaussian  distribution 

C  function  for  a  zero  mean,  unit  variance  random  variable. 

C 

C  ROUT I  NEC S)  ACCESSED  OR  CALLED  BY  THIS  ROUTINE 

C 

C  NONE 

C 

C  ARGUMENT(S)  REQUIRED  FROM  THE  CALLING  ROUTINE 

C 

C  X  value  at  which  to  evaluate  the  error  function 

C 

C  ARGUMENT(S)  SUPPLIED  TO  THE  CALLING  ROUTINE 

C 

C  NONE 

C 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NOTE  1:  This  routine  us^s'  two  different  approximations 
to  compute  Q(x)  depending  on  the  value  of  x. 

For  0  <  X  <  2.326  ,  see  Note  2. 

For  X  >  2.326  ,  see  Note  3. 

NOTE  2:  For  values  of  X  in  the  range  0  <  X  <  2.326, 

(I.e.,  for  .5  >  Q(x)  >  .01  ),  equation  26.2.17 
on  page  932  of  reference  1  is  used  to  evaluate  Q(x). 

NOTE  3:  For  values  of  X  >  2.326  (i.e.,  for  Q(x)  <  .01  ) 
equation  13  on  page  641  of  reference  2  is  used 
to  eval uate  Q(x) . 

NOTE  4:  These  approximations  insure  that  the  error  in  the 
value  of  Q(x)  Is  less  than  0.001  %  . 

NOTE  5:  This  routine  works  for  positive  and  negative  values 
of  X.  Note  that  Q(-x)  =  1  -  Q(x). 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


REFERENCE  1:  Abramowitz  and  Stegun,  HANDBOOK  OF  MATHEMATICAL 

FUNCTIONS,  National  Bureau  of  Standards,  1964. 

REFERENCE  2:  P.  0.  Borjesson,  C.  E.  Sundberg,  SIMPLE 

APPROXIMATIONS  OF  THE  ERROR  FUNCTION 'Q(x)  FOR 
COMMUNICATION  APPLICATIONS,  IEEE  Transactions 
on  Communications,  March  1979,  pp  639-643. 


c********************************************«*****»*********»******#**» 

C 

REAL  FUNCTION  Q  (X) 

C 

DOUBLE  PRECISION  F(5),  P,  T,  ZEXP,  QTEMP,  SQT2PI 
REAL  A,  B,  X,  XPOS,  DEN 
C 

DATA  SQT2PI  /2. 506628275/,  P  /0.2316419/ 

DATA  F  /O. 319381530,  -.356563782,  1.781477937,  -1.821255978, 

&  1 .330274429/ 

DATA  A  /0.279/,  B  /7.123/ 

C 

C********************************* 

C 

C  Initialize  Variables 

C 

C 


QTEMP  =0.0 
XPOS  =  ABS(X) 

ZEXP  =  EXP(-(X**2)/2.)/SQT2PI 
C 

c 

C  Determine  range  of  X,  and  compute  Q(x) 

C  using  the  appropriate  approximation 

C 

C 

IF  (XPOS.LE.2.326)  THEN 

T  =  1  .  /  ( 1  .  +  P  *  XPOS ) 

DO  I  =  1,  5 

QTEMP  =  QTEMP'  +  F(l)  *  (T**l) 

END  DO 

QTEMP  =  QTEMP  *  ZEXP 

ELSE 

DEN  =  (1.  -  A)  *  XPOS  +  A  *  SQRT(B  +  XP0S**2) 
QTEMP  =  ZEXP  /  DEN 

END  IF 
C 

C****************************************** 

c 

C  Check  for  negative  values  of  X,  and 
C  adjust  Q(x)  If  necessary 
C 
C 


IF  (X.GE.O.)  THEN 
Q  =  QTEMP 

ELSE 

Q  =  1 .0  -  QTEMP 
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END  IF 

RETURN 

END 
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APPENDIX  D 


DERIVATION  OF  EXPRESSION  FOR  Pg(T,  t+a) 


The  derivation  of  the  expression  for  Pg(T,  t+a),  equation  (3-30), 
is  given  in  [29]  and  is  included  in  this  appendix  for  completeness.  Also 
considerations  leading  to  the  bound  for  Pg(A)  in  (3-31)  are  given. 

Recall  Pg(T,  t+a)  represents  the  minimum  attainable  probability 
of  error  for  the  binary  decision  problem  of  deciding  which  of  the  two 
hypothesized  delay  values,  x  or  t+a,  is  correct.  For  large  BT, 

Pg(T,  t+a)  can  be  approximated  by  the  Chernoff  formula  [23,  p,  125, 


eq,  (484)] 


Pe("’  ■  =iexp[l(q^)  +  q^  l"/2]  Q[qg  n/tI^] 

+  \  exp[l(q^)  +  (1  -  q^)2  l"(qg)/2]  Q[(l  -  q^)  v/l"(qg)]  (D-la) 


where 


oo 


(D-lb) 


and  where  q^  is  the  value  for  which  l'(q)  =0  (note,  l'(q)  =  dl/dq  and 


p  p 

l"(q)  =  d'^l/dq  ),  R  is  the  observation  vector  made  up  of  the  Fourier 


coefficients  of  the  received  signals,  r^{t)  and  r2(t),  and  p(^|t)  is 
the  conditional  pdf  of  R  given,  as  discussed  in  Appendix  A, 
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Substituting  the  expressions  for  p(R|t)  and  P(^|t+a)  from  (A-6)  into 
(D-lb)  and  carrying  out  the  integration  yields 


l(q)  =  -  jE 


q  In  [detEi^i^  ]  +  (1-q)  In  CdetEk|^+^] 


+  In  [det  (q^klT  ^  Skk+A 


(D-2) 


As  noted  in  (A-8),  det^i^  |  ^  is  independent  of  t.  Recall  from  (A-7)  that 


ss'  o'  nn'  o' 


-jkw„T 


aG^s(kWo)  e 


jkw^x 


aGss(kWo)  e 


o  G  (kw  )  +  G  (kw  ) 
ss  o'  nn'  O' 


(D-3a) 


thus 


-1 


2k  lx 


■  ? 


.-a  G3^(kw^)  e 


-jkw  X 


(D-Sb) 


Substituting  (D-3)  into  (D-2)  and  simplifying  yields 


l(q)  =  -  7  E  In  [1  +  4q(l-q)  Y(kf^,A)] 
^  k  ° 


(D-4) 


where  Y(f,A)  is  given  in  (3-30e).  Then  taking  the  first  and  second 
derivatives  of  l(q)  with  respect  to  q  gives 


^2  (l-2q)  Y(kf^,A) 

2  1  +  4q(l-q)  Y(kfQ,A) 


(D-5a) 


and 
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(D-5b) 


l"(q)  = 


k 


4  Y(kfQ,A) 

1  +  4q(l-q)  Ylkf^^A) 


8(l-2q)^  Y^(kfQ,A) 

[1  +  4q(l-q)  T(kf^,A)]2. 


Note  that  l'(qQ)  =  0  when  q^  =  1/2.  For  large  BT,  the  summations  in 
(D-4)  and  (D-5)  can  be  approximated  as  integrals  using  the  large  T 
approximation  as  in  Appendix  A.  Evaluating  l(q)  and  r'(q)  at  q^  =  1/2 
yields 

CP 

1(1/2)  =  -  J  ^  Y(f,A)]  df 

—  oo 

and 


l"(l/2)  =41^^ 


— iLL;A) 

1  +  Y(f,A) 


df 


Substituting  q^  =  1/2  and  (D-6)  into  (D-la)  and  simplifying  yields  the 
expression  for  Pg(T,  t+a)  given  in  (3-30),  namely 


Pg(  T,  t+a)  =  Pg(A) 

where 

Pg(A)  =  exp(a(A)  +  b(A))  Q('sy2b(A)), 

and  where 

oo 

a(A)  =  -  T y*  In  (1  +  Y(f,A))  df  , 
0 


b(A)  =1  f  Y(f,A)/(l  +  Y(f,A))  df  , 
0 


and 


Y(f»A)  = 


G  (f)^  sin^irfA 
ss 


IG 


"l"l 


"l"l 


(f)  G„  „  (f) 


n  2  ^2 


(D-6a) 

(D-6b) 

(D-7a) 

(D-7b) 

(D-7c) 

(D-7d) 

(D-7e) 
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Now  to  obtain  the  bound  for  Pg(A)  given  in  (3-31),  let  y  =  Y(i'»A)  and 
note  that  0  <  y/(1''’y)  <  1*  Also,  note  that 


1n(l+Y) 


thus. 


From  (D-8b)  it  is  easily  seen  that 


0  >  a(A)  +  b(A) 

nn(l  +  Y(f,A)) 


oo 

.-t/[ 

0  I- 


Y(f,A)  ~| 

■  i'  +  Y(f,A)J 


df 


>  —  a)  dF  . 

~  "^0 


(D-8a) 


(D-8b) 


(D-8c) 


(D-9) 
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Also,  since  y/(1'^y)  <  y\ 


on 

0  £  b(A)  1  ^  y*  Y(f»A)  df  . 

0 

Substituting  (D-9)  and  (D-10)  into  (D-7b),  immediately  yields  (3-31), 


“■  00  ^ 

00  — 

/  r  \ 

Pg(A)  >  exp 

4^  df 

•  Q 

\21  J  Y(f,A) 

(D-10) 


(D-11) 


120 


APPENDIX  E 


APPROXIMATION  FOR  THE  THRESHOLD  SNR 


The  approximations  leading  to  the  expression  for  the  threshold  SNR  of 
(4-2)  are  given  in  this  appendix.  The  initial  assumption  is  that  the 
threshold  SNR  can  be  closely  approximated  by  the  SNR  which  satisfies  the 
condition 

2  _  9  2 
‘^ZZLB  “  ^‘^CRLB  • 

From  (3-61) 

[(nz)^  -  1]  Q(nz)  -  +  1 


2  1 

®ZZLB  “  7T 
2n 


+  2Pg  D^/S 


(E-2a) 


where 


nz 


«V1  +  SNR'  -  1 
^/l  +  SNR' 


1/2 


(E-2b) 


and 


P 


e 


exp 


-BT 


?  1  n  f 

1  +  VI  +  SNR'  \ 

C  1  1 

ro 

VI  +  SNR'  -  1 
Vl  +  SNR' 


Q(nz) . 


For  large  BT,  the  threshold  SNR  is  small,  then  assume  that  SNR'  <<  1. 
assumption  allows  the  following  approximations  to  be  made. 


(E-2c) 

This 
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(E-3a) 


1  +  n/1  +  SNR' 

1  +  n/1  +  SNR' 

1  +  n/T 

+  SNR' 

2 

2 

_1  +  n/T 

+  SNR' 

=  1  + 


SNR' 


2(1  +  sll  +  SNR') 


=  1 


+ 


SNR 

4 


and  similarly, 

n/1  +  SNR'  -  1  ^  N/f  +  SNR '  -  1  n/1  +  SNR'  +1 
n/1  +  SNR'  n/1  +  SNR'  [Nyn^NR'  +  1  _ 

^  _ SNRJ _ 

1  +  SNR'  +  ^yl  +  SNR' 

SNR' 

"  ~r~ 

Substituting  (E-3a)  and  (E-3b)  -into  (E-2a)  and  (E-2b)  yields 


nz  =  n/BT  SNR' 
and 

Pg  =  exp  |-BT  [2  1n(l  +  SNR'/4)  -  SNR'/2]^  Q(^/BT  SNR'). 


Note  that  1n(l  +  x)  =  x  for  small  x,  thus 


Pg  =  Q(n/BT  SNR'). 


Further  suppose  that  BT  SNR'  is  large  enough  so  that 


e-(nz)^/2  ^  g-(BT  SNR')/2  ^  q 


(E-3b) 

(E-4a) 

(E-4b) 

(E-5) 

(E-6) 
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Substitution  of  (E-4)-(E-6)  into  (E-2a)  yields 


a^ZZLB  “T  [(BT  SNR'  -  1)  Q(n/BT  SNR'  )  + 

2n 

+  2  QWBT  SNR')  D^/3  .  (E-7) 

Recall  from  (3-63)  and  (3-64)  that 

‘'CRLB  =  TT 
4n 

where 

=  2ii^  SNR'  B^T/3  .  (E-8b) 

Substituting  (E-7)  and  (E-8a)  into  (E-1)  and  simplifying  yields 

(BT  SNR'  -  1)  Q(  n/BT  SNR')  +  4Q(  n/BT  SNR')  Dq  ti^/3  =  |  (E-9a) 

and  substituting  for  n  from  (E-8b)  gives 

BT  SNR'  Q(n/BT  SNR'  )  [1  -  1/(BT  -  SNR')  +  8ir^B^D^/9]  = .  '  (E-9b) 

In  the  derivation  of  the  ZZLB  it  was  assumed  that  >>  1/2B,  thus, 

BD^  >>  1  and  (E-9)  becomes 


BT  SNR'  Q(n/BT  SNR')~=  - 

16 

0 


(E-10) 


Finally  solving  (E-10)  for  SNR'  yields  the  approximate  expression  for 


SNR'^^  of  (4-2), 


SNR'th 


(E-lla) 
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where 


F(x)  =  Q(x), 


(E-llb) 


and  F;^(- 


denotes  the  larger  of  the  two  solutions  of  the  inverse  of  F(*)- 
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APPENDIX  F 


THE  ML  ESTIMATOR  -  TIME  VARYING  TIME  DELAY 


This  appendix  briefly  summarizes  the  derivation  qiven  by  Knapp  and 
Carter  in  [46]  for  the  ML  estimator  for  the  time-varying  TDE  problem  of 
(5-1).  The  procedure  is  very  similar  to  that  used  to  derive  the  ML 
estimator  for  the  case  of  a  fixed  time  delay  as  discussed  in  Appendix  A, 
and  similar  notation  is  used  here.  The  ML  estimate  for  the  parameters 

62  and  D  is  obtained  by  maximizing  the  conditional  pdf  p(^  |  bj^,  b2,  t  ), 
where  b^,  b2,  and  t  represent  hypothesized  values  for  Bp  B2»  and 
D,  respectively.  The  conditional  pdf  can  again  be  expressed  as 


p(R  I  b^,  )  =  c  exp(-  ^  J) 


{F-2a) 


where 


k=-N 


and 


N 


-1 


{F-2c) 


However,  for  the  signal  of  (5-1) 


=  E  CR(k)  R*'(k)  b^,  b2,T  ] 
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-jkw^bx 

(bkWg) 


=  y  Q,  » 


(F-3) 


where  b  =  b2  /b^. 


In  practice,  the  power  spectral  matrix  £  must  be  estimated  from  the 
data.  Since  b^,  b2  =  1,  ^  may  be  replaced  by 


Q(f)  =  Q(f) 


Gr  p  (n 

'^ri 


G*  ^  (f) 


Gr  p  (n 


G  (f) 
'“3'^3 


(F-4a) 


where 


r3(t)  =  r2(t/b)/  b 


(F-4b) 


and  hence 


'“l'“3 


(f) 


g-j2irfbT 


(F-4c) 


That  is,  '"3(t)  represents  a  time  scaled  version  of  r2(t)  where  the  time 
scaling  compensates  for  the  relative  time  compression  between  and 
r2(t).  Then,  the  quantity  c  in  (F-2)  is  independent  of  x  and  is 


126 


approximately  independent  of  and  b2  (since  b2  =  1),  so  that 
maximizing  p(£  |  b^^,  b2,  t)  is  equivalent  to  minimizing  J  in  (F-2b). 

For  large  T, 


J 


R*'(f)  Q(f)‘^  R(f) 


df . 


(F-3) 


By  substituting  for  2(f)~^  and  manipulating  (F-3),  it  can  be  shown  that 
minimizing  J  is  equivalent  to  maximizing 


00 

R^(f)  R2(6f)  W(f)  df  ,  (F-4a) 

—  CO 

or  equivalently. 


00 


—  00 


Rl(f)  R3(f)  W(f)  df 


(F-4b) 


where 


W(f)  = 


Cr  r 


Gr  r 


(1  -  r  (f)) 


(F-4c) 


As  noted  in  Chapter  5,  the  ML  estimates  of  6  =  ^2^^!  ^  correspond 

to  the  values  of  b  and  t  for  which  J'  (or  J")  is  a  maximum. 
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APPENDIX  G 


SIMULATION  DETAILS  -  TIME-VARYING  TIME  DELAY 

A  computer  simulation  routine  implementing  the  compensation  technique 
described  in  section  5.3  has  been  written  and  preliminary  simulation  results 
have  been  obtained  as  discussed  in  section  5.4.  The  signal  and  noise 
sequences  were  generated  in  the  same  manner  as  described  in  Appendix  B. 
However,  in  this  case,  a  time  shifted  version  of  the  signal  sequence  was 
generated  using  the  technique  proposed  by  Youn  [14,  56].  Program  listings 
for  the  simulation  main  routine,  DOPPLER,  and  the  subroutine  SHIFT,  which 
generates  the  time  shifted  signal  are  included  in  this  appendix.  The 
portions  of  the  main  routine  which  specifically  relate  to  the  compensation 
technique  have  been  boxed.  Otherwise  the  processing  is  essentially 
identical  to  that  for  the  case  of  a  fixed  time  delay.  The  functions  and 
subroutines  called  by  the  main  routine  are  given  in  the  header  of  the 
program  listing,  however,  the  subroutine  listings  (other  than  SHIFT)  have 
not  been  included.  Variables  are  also  documented  in  the  listings.  However, 
it  is  worth  noting  that  a  maximum  value  on  the  delay  value  (MAXDEL)  is  set 
in  the  program.  When  the  time  delay  exceeds  this  value  in  absolute  value, 
the  sign  of  time  delay  rate  is  changed  (D_RATE=-D_RATE) . 
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****** ************************************************************* 

c 

C  DOPPLER  COMPENSATION  SIMULATION  ROUTINE 

C 

C  VAX-11  FORTRAN  SOURCE  FILENAME:  DOPPLER. FOR 

C 

C  DEPARTMENT  OF  ELECTRICAL  ENGINEERING  KANSAS  STATE  UNIVERSITY 

C 

C  REVISION  DATE  PROGRAMMER(S) 

C  -  -  - 

C  01.0  JULY  1983  KENT  SCARBROUGH 

C 

Q*********************************************************************** 

c 

C  PURPOSE 

C 

C  This  routine  implements  a  simulation  to  test 

C  a  frequency  domain  technique  to  compensate  for 

C  the  "doppler"  shift  which  occurs  due  to  relative 

C  motion  between  source  and  sensors  in  the  time 

C  delay  estimation  problem. 

C 

C  ROUTINE(S)  CALLED  BY  THIS  ROUTINE 

C 

C  BUTTER,  FFTCX,  F0RD2,  GCC,  GWNSEQ,  HAMMING,  HEADER, 

C  PEAK,  SGOPEN,  SGTRAN,  REPLY,  RSEED,  RVINIT,  SHIFT, 

C  SSEED,  WAML,  WMSC,  WOPT,  WPHT,  WSCC,  WSCT 

C 

C  OTHER  ROUT I  NEC S)  ACCESSED  BY  TH I S  ROUTINE 

C 

C  FFS,  FFT842,  LPDES,  HPDES,  BPDES,  BSDES,  COEFF,  GAUSS 

C 


Q*********************************************************************** 

c*********************************************************************** 

c 

c 


REAL  R1 (32*1024),  R2(32*1024),  SGI (32*1 024+61 ) ,  SG2(32*1024) 
REAL  R1W(75),  R2W(75),  SG1W(75),  SG2W(75),  PAST(75) 

REAL  WORKX(2050),  WORKY(2050) 

REAL  RXX(2048),  RYY(2048),  Wl NDOW( 1024) ,  WEI GHT( 1025)  - 
REAL  GXX(1025),  GYY(1025),  GXYR(1025),  GXYI(1025) 

REAL  GXYRSVC 1 025,21 ),  GXYI SVC  1 025,21 ) 

REAL  TDELAY(500),  TDRATE(500) 

REAL  SCCTAUC  500),  SCTTAUC  500),  AMLTAUC  500) 

REAL  MSCTAUC  500),  OPTTAUC  500),  PHTTAUC  500) 

REAL  SCCTDTC  500),  SCTTDTC  500),  AMLTDTC  500) 

REAL  MSCTDTC  500),  OPTTDTC  500),  PHTTDTC  500) 

REAL  SCCPK,  SCTPK,  AMLPK,  MSCPK,  OPTPK,  PHTPK 

REAL  DELDOT,  DEL_F,  DEL_T,  D_RATE,  DM  IN,  FNORM,  FS 
REAL  NSFR(?L,  NSFR(?U,  SGFR(?L,  SGFR(?U,  PKVAL,  SCALE 
REAL  SNR1,  SNR2,  SSTP,  TIME,  TDEL1 

INTEGER  ADELAY,  MAXDEL,  FSECT,  NSECT,  IPNT,  LPNT 
INTEGER  PKLOC,  NCMP,  NTRL,  NPTS,  NPTP1 ,  NPTM2,  NPTD2 
INTEGER  NFFT,  NSEG,  NS  IN,  TSEG,  TPTS,  WIDTH,  UN  I  TP 
INTEGER  NSEED,  SEED(38) 

LOGICAL  see,  SCOT,  AML,  MSC,  OPT,  PHAT,  RDSEED,  SVSEED 
LOGICAL  FILTER,  RESET,  NOSET 
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CHARACTER  FTYPE*2,  SNRDB*3,  SNRTMP*2,  FCTN»40,  TNAM*40 
C 
C 

C 

C  INITIALIZATION 
C 


DATA 

TDEL1 

/  0/, 

MAXDEL 

/  15/, 

DELDOT  / 

2.5E-4/ 

DATA 

D  RATE 

/  2.5E-4/, 

FS 

/  2048./, 

FTYPE 

/ 

'LP*/ 

DATA 

FSECT 

/  21, 

NSFRQL 

/  0.0/, 

NSFRQU  / 

100./ 

DATA 

NCMP 

/  21, 

NTRL 

/  51/, 

NPTS 

/ 

512/ 

DATA 

SGFRQL 

/  0.0/, 

SGFRQU 

/  100./, 

SNR1 

/ 

-3.0/ 

DATA 

SNR2 

/  -4.0/, 

SSTP 

/  -3.0/, 

TIME 

/ 

4.0/ 

DATA 

NWT 

/  51/, 

WIDTH 

/  256/ 

DATA 

see 

/  .TRUE./, 

SCOT 

/  .TRUE./, 

AML 

/. 

FALSE./ 

DATA 

MSC 

/.FALSE./, 

OPT 

/.FALSE./, 

PHAT 

/. 

FALSE./ 

DATA 

F 1 LTER 

/  .TRUE./, 

RESET 

/  .TRUE./, 

NOSET 

/. 

FALSE./ 

DATA 

RDSEED 

/.FALSE./, 

SVSEED 

/.FALSE./, 

NSEED 

/ 

38/ 

C 

******************** ************* ************** 

c 

C  FORMAT  STATEMENTS 

C 

1  FORMAT(A) 

3  FORMAT(A, I10,A) 

4  F0RMAT(A,F10.2,A) 

5  F0RMAT(A,8X,A) 

7  F0RMAT(I2.2) 

C 

Q*** *************************************************** ******* 

C 

C  CALCULATE  CONSTANTS  -  INITIALIZE  PROGRAM  VARIABLES 

C 

FNYQ  =  FS  /  2. 

NPTP1  =  NPTS  +  1 
NPTM2  =  NPTS  *  2 
NPTD2  =  NPTS  /  2 
NFFT  =  NPTM2 
TPTS  =  IFIX(TIME  *  FS) 

NSEG  =  TPTS  /  NPTS 
TSEG  =  2  *  NSEG  -  1 
C 

DEL_F  =  FNYQ  /  FLOAT (NPTS) 

DEL_T  =  TIME  /  FLOAT(TSEG) 

TWOPI  =  2.0  *  3.141593 
NCMP21  =  NCMP  *2+1 
NWT2  =  ( NWT  -  1  )  /  2 
C 

IPNT  =  NPTS  -  WIDTH/2  +  1 
LPNT  =  NPTS  +  WIDTH/2 
DMIN  =  FLOAT( IPNT  -  NPTS) 

C 

IF  ((FTYPE.EQ. 'LP').OR.(FTYPE.EQ. 'HP'))  THEN 
NSECT  =  2*FSECT 
ELSE  IF  (FTYPE.EQ. 'BP' )  THEN 
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FSECT  =  FSECT  /  2 
NSECT  =  4*FSECT 

ELSE 

FILTER  =  .FALSE. 

END  IF 
NSAV  =  75 

IF  ((NSFRQU.EQ.SGFRQU).AND.(NSFRQL.EQ.SGFRQL))  RESET  =  .FALSE. 

COMPUTE  HAMMING  WINDOW 

CALL  HAMM  I  MG ( W I NDOW , NPTS , WPOWER ) 

FNORM  =  TSEG  *  NPTS  *  WPOWER 

INITALIZE  RANDOM  NUMBER  GENERATOR 

IF  (RDSEED)  THEN 

CALL  SGOPEN ( 98 , ' READ ' , ' NOPROMPT ' , ' SEED . GCC ' , ' I NTEGER ' , NSEED ) 
CALL  SGTRAN(98,'READ',' INTEGER' , SEED, NSEED) 

CALL  RSEED(SEED(1 ),SEED(7)) 

END  IF 

UN  I  TP  =  99 


************************************************************* 


ECHO  CHECK  OF  PARAMETERS 


CALL  HEADER(UNITP 

, '***  MOVING  SOURCE  SIMULATION  *** ' 

) 

WRITE(UNITP,3)  * 

Number  of  Trials 

f 

f 

NTRL 

WRITE(UNITP,4)  ' 

Sampling  Frequency 

f 

f 

FS  , 

f 

Hz  ' 

WRITE(UNITP,4)  ' 

Observation  Time 

f 

f 

TIME  , 

f 

seconds ' 

WRITE(UNITP,3)  ' 

Segment  Length 

f 

f 

NPTS  , 

f 

points  ' 

WRITE(UNITP,3)  ' 

FFT  Size 

f 

f 

NFFT  , 

f 

points  ' 

WRITE(UNITP,3)  ' 

Correlation  Window 

1 

P 

WIDTH, 

f 

points  ' 

WRITECUNITP,*) 

WRITE(UNITP,5)  ' 

Fi 1  ter  Type 

f 

P 

FTYPE 

WRITE(UNITP,3)  ' 

Fi  1  ter-  Order 

f 

P 

NSECT 

WRITE(UNITP,4)  ' 

Lower  Cutoff  (noise) 

f 

P 

NSFRQL 

P 

'  Hz' 

WRITE(UNITP,4)  ' 

Upper  Cutoff  (noise) 

f 

P 

NSFRQU 

P 

'  Hz' 

WRITE(UNITP,4)  ' 

Lower  Cutoff  (signal) 

f 

P 

SGFRQL 

P 

'  Hz' 

WRITE(UNITP,4)  ' 

Upper  Cutoff  (signal) 

f 

P 

SGFRQU 

P 

'  Hz' 

WRITE(UNITP,4)  ' 

First  SNR  Value(dB) 

1 

P 

SNR1 

P 

'  dB' 

WRITE(UNITP,4)  ' 

Last  SNR  Value(dB) 

1 

P 

SNR2 

P 

'  dB' 

WRITE(UNITP,4)  ' 

SNR  Steps ize(dB) 

1 

P 

SSTP 

P 

'  dB' 

TEMPI  =  D_RATE  * 

1 .0E4 

TEMP2  =  DELDOT  * 

1 .0E4 

WRITECUNITP,*) 

WRITE(UNITP,3)  ' 

Number  of  Comp. 

1 

P 

NCMP 

WRITE(UNITP,4)  ' 

Initial  Delay  Value 

P 

TDELl 

WRITE(UNITP,3)  ' 

Maximum  Delay  Value 

} 

P 

MAXDEL 

WRITE(UNITP,4)  ' 

Initial  Delay  Rate 

P 

TEMPI , 

E-4' 

WRITE(UNITP,4)  ' 

Assumed  Rate  Steps ize 

1 

P 

TEMP2, 

» 

E-4' 

WRITE(UNITP,3)  ' 

No.  of  Filter  Coeffs 

\ 

P 

NWT 
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c 


WRITE(UNITP,*) 

IF  (RDSEED)  WRITE(UNITP,*) 

&  'Initial  Seed  Values  Read  From  File:  SEED.GCC 
IF  (SVSEED)  WRITE(UNITP,*) 

&  'Final  Seed  Values  Output  To  File:  SEED.GCC 
C 

WRITE(UNITP,*) 

C 

Q*»*»******************************************************************* 

Q**»******************************************************************** 

C 

C  MAIN  LOOP  :  EXECUTED  ONCE  FOR  EACH  SNR  VALUE 
C 

DO  SNR  =  SNR1 ,  SNR2,  SSTP 
C 

ASNR  =  10.**(SNR/10. ) 

SCALE  =  1.  /  SORT(ASNR) 

C 

ISNR  =  IFIX(SNR) 

ITMP  =  ABS(  ISNR) 

WRITE(SNRTMP,7)  I  TOP 
IF  ( ISNR.LT.O)  THEN 

SNRDB  =  'N'  //  SNRTMP 

ELSE 

SNRDB  =  'O'  //  SNRTMP 
END  IF 
C 

Q************************************************************* 

c************************************************************* 

c 

C  SECONDARY  LOOP  :  EXECUTED  ONCE  FOR  EACH  TOIAL 
C 

DO  ITRIAL  =  1,  NTRL 
C 
C 

Q*************************************************** 

Q*************************************************** 

C 

C  GENERATE  RECEIVED  S I GNAL' SEQUENCES 

C 

C  GENERATE  WH ITE,GAUSS I  AN  NOISE  SEQUENCES 

C 

CALL  GWNSEQ(R1 ,TPTS, 0.0,1 .0) 

CALL  GWNSEQ(R2,TPTS,0.0,1 .0) 

IF  ( ITRIAL.EQ.1 ) 

&  CALL  GWNSEQ(PAST,NWT2, 0.0,1  .0) 

CALL  GWNSEQ(SG1 ,TPTS, 0.0,1 .0) 

C 

C*************************************************** 

c 

C  FILTER  SEQUENCES  TO  OBTAIN  DESIRED  SPECTRA 

C 

C 

IF  (FILTER)  THEN 

CALL  BUTTER(R1 ,TPTS,R1W,NSFRQL,NSFRQU,FS, 

&  FTYPE,FSECT, RESET) 
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CALL  BUTTER ( R2 , TPTS , R2W, NSFRQL , NSFRQU, FS , 

&  FTYPE,FSECT, NOSET) 

IF  ( ITRIAL.EQ.I )  THEN 

CALL  BUTTER ( PAST , NWT2 , SGI W, SGFRQL , SGFRQU , FS , 

&  FTYPE,FSECT, RESET) 

CALL  BUTTERCSGI , TPTS, SGI W, SGFRQL, SGFRQU, FS, 

&  FTYPE,FSECT, NOSET) 

ELSE 

CALL  BUTTERCSGI , TPTS, SGIW, SGFRQL, SGFRQU, FS, 

&  FTYPE,FSECT, RESET) 

END  IF 
END  IF 

C 

C  SHIFT  SIGNAL  SEQUENCE,  APPEND  PAST  VALUES 

C  AND  SAVE  ’NEW  PAST  VLUES 

C 

DO  IPT  =  TPTS,  1,  -1 
JPT  =  IPT  +  NWT2 
SGI ( JPT)  =  SGI (IPT) 

END  DO 
C 

DO  IPT  =  1,  NWT2 

JPT  =  TPTS  +  IPT 
SGI (IPT)  =  PAST( IPT) 

PAST( IPT)  =  SGI (JPT) 

END  DO 


c 

c 

c 

c 

c 


C 


SCALE  NOISE  SEQ  TO  OBTAIN  DESIRED  SNR, 
ADD  SIGNAL  TO  NOISE  SEQUENCES 


DO  IPT  =  1, 
R1 ( IPT) 
R2( IPT) 
END  DO 


TPTS 

=  SCALE  *  R1 ( IPT)  +  SGI ( IPT) 
=  SCALE  *  R2( IPT)  +  SG2( IPT) 


c 

c 

c 

c 

c 


IMPLEMENT  GCC  PROCESSOR 

ZERO  SUMMING  ARRAYS  &  WORK  ARRAYS 
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CALL  RVINIT(GXX,NPTP1 ,0.0) 

CALL  RVINIT(GYY,NPTP1 ,0.0) 

CALL  RVINIT(GXYR,NPTP1 ,0.0) 

CALL  RVINIKGXYI  ,NPTP1 ,0.0) 

CALL  RVINIKGXYRSV,  1025  »  21,0.0) 

CALL  RVINIT(GXYISV,1025  *  21,0.0) 

CALL  RVINIT(W0RKX,NPTM2+2,0.0) 

CALL  RVINlT(W0RKY,NPTM2+2,0.0) 

C 

SCCPK  =0.0 
SCTPK  =  0.0 
AMLPK  =  0.0 
MSCPK  =0.0 
OPTPK  =  0.0 
PHTPK  =0.0 

Q ********* ****************************************** 

C 

C  LOOP  EXECUTED  ONCE  FOR  EACH  DATA  SEGMENT 

C 

DO  ISEG  =  1,  TSEG 
C 

0*************************************************** 

C 

C  WEIGHT  DATA  ARRAYS  WITH  HAMMING  WINDOW 

C  AND  PAD  WITH  ZEROES  (DOUBLE  LENGTH) 

C 

DO  IPT  =  1,  NPTS 

IND1  =  (ISEG  -  1 )  *  NPTD2  +  IPT 
RXX(IPT)  =  RKINDI)  *  WINDOW(IPT) 

RYY(IPT)  =  R2(IND1)  *  WINDOW(IPT) 

IND2  =  NPTS  +  IPT 
RXX(IND2)  =  0.0 
RYY(IND2)  =  0.0 
END  DO 
C 

0*************************************************** 

C 

C  COMPUTE  FORWARD  FFT  AND  AVERAGE  THE 

C  PERIODOGRAM  AUTO-SPECTRAL  ESTIMATES 

C  ■ 

CALL  FFTCX(RXX,RYY,NFFT, 'FORWARD' , 'NONORM' ) 

CALL  F0RD2 ( RXX, RYY , NFFT, WORKX , WORKY ) 

C 

GXXd)  =  GXXd)  +  RXXd)  »  RXXd) 

GYYd)  =  GYYd)  +  RYYd)  *  RYYd) 

GXYRd)  =  GXYRd)  +  RXXd)  *  RYY(I) 

DO  IPT  =  2,  NPTS 

JPT  =  2  »  ( IPT-1 ) 

JP1  =  JPT  +  1 

GXXdPT)  =  GXXdPT)  +  RXX(JPT)**2  +  RXX(JP1)**2 
GYYdPT)  =  GYYdPT)  +  RYY(JPT)**2  +  RYY(JP1)**2 
GXYRdPT)  =  RXX(JPT)  *  RYY(JPT) 

&  +  RXX( JP1 )  *  RYY(JP1 ) 

GXYldPT)  =  RXX(JPT)  *  RYY(JPI) 

&  -  RXX( JP1 )  »  RYY(JPT) 

END  DO 

GXX(NPTPI)  =  GXX(NPTPI)  +  RXX(NPTM2)  *  RXX(NPTM2) 
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GYY(NPTPI)  =  GYY(NPTPI)  +  RYY(NPTM2)  *  RYY(NPTM2) 
GXYR(NPTPI)  =  GXYR(NPTPI)  +  RXX(NPTM2)  *  RYY(NPTM2) 


^ - 

Q*************************************************** 

C 

C  COMPENSATE  FOR  TIME  SHIFT  BETWEEN  DATA  SEGMENTS 

C  &  AVERAGE  COMPENSATED  CROSS-SPECTRAL  ESTIMATES 

C 

DO  JCMP  =  1 ,  NCMP21 

ICMP  =  JCMP  -  (NCMP  +  1 ) 

TAUDOT  =  FLOAT(ICMP)  *  DELDOT 

RSEG  =  FLOAT(ISEG)  -  .5 

ARGTTIP  =  TWOPI*  RSEG  *  DEL_T  *  TAUDOT 

GXYRSV( 1 ,JCMP)  =  GXYRSV( 1 , JCMP)  +  GXYR(I) 

GXYISV(1,JCMP)  =  0.0 

DO  IPT  =  2,  NPTS 

FREQ  =  FLOATdPT  -  1)  *  DEL_F 
ARG  =  FREQ  *  ARGTOP 
JPT  =  2  *  (IPT-1) 

JP1  =  JPT  +  1 

GXYRSV( IPT, JCMP)  =  GXYRSV ( I PT, JCMP)  + 

&  GXYR(IPT)  *  COS(ARG)  -  GXYKIPT)  *  SIN(ARG) 

GXYISV( IPT, JCMP)  =  GXYISV( IPT, JCMP)  + 

&  GXYR(IPT)  *  SIN(ARG)  +  GXYKIPT)  *  COS(ARG) 

END  DO 

ARG  =  FNYQ  *  ARGTTIP 

GXYRSV ( NPTP 1 , JCMP)  =  GXYRSV ( NPTP1 , JCMP)  + 

&  GXYR(NPTPI)  *  COS(ARG) 

GXYISV(NPTP1 , JCMP)  =  0.0 
END  DO 

-G - - 

Q^i*^i*************************:*:**^(:******************* 


C 

C  GO  BACK  FOR  NEXT  DATA  SEGMENT 

C 

END  DO 


C 


Q*************************************************** 


c 

c 

c 

c 


c 


LOOP  EXECUTED  ONCE  FOR  EACH  COMPENSATION  VALUE 


DO  JCMP  =  1 ,  NCMP21 

ICMP  =  JCMP  -  (NCMP  +  1 ) 

TAUDOT  =  FLOAT(ICMP)  *  DELDOT 

DO  IPT  =  1,  NPTP1 

GXYR(IPT)  =  GXYRSV (IPT, JCMP) 
GXYKIPT)  =  GXYISV(  IPT,JCMP) 
END  DO 


C 

C  NORMALIZE  SPECTRAL  ESTIMATES  TO  ACCOUNT  FOR: 

C  1)  POWER  REDUCTION  DUE  TO  WINDOWING 

C  2)  SUMMING  TSEG  ESTIMATES 
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3)  NUNBER  OF  POINTS/SEGMENT 
NOTE:  NO  NORMALIZATION  IS  DONE  IN  FFT'S 

IF  (JCMP  .EQ.  1)  “mEN 
DO  IPT  =  1,  NPTP1 

GXX(IPT)  =  GXX(IPT)  /  FNORM 
GYY(IPT)  =  GYY(IPT)  /  FNORM 
END  DO 
END  IF 

DO  IPT  =  1 ,  NPTP1 

GXYR(IPT)  =  GXYR(IPT)  /  FNORM 
GXYKIPT)  =  GXYKIPT)  /  FNORM 
END  DO 


COMPUTE  GCC  FUNCTIONS  FOR  EACH  COMPENSATION  VALUE 
DETERMINE  GLOBAL  PEAK  IN  TAU-TAUDOT  COORDINATES 


IF  (SCO  “mEN 

CALL  RVINIT(WEIGHT,NPTP1 ,1 .0) 

CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 
CALL  PEAK(WORKX, IPNT,LPNT,PKVAL,PKLOC) 

IF  (SCCPK  .LT.  PKVAL)  “mEN 
SCCPK  =  PKVAL 
IND  =  PKLOC  -  IPNT  +  1 
SCCTAUC ITRIAL)  =  FLOATCIND  -  WIDTH/Z) 
SCCTDTC ITRIAL)  =  TAUDOT 
END  IF 
END  IF 

IF  (SCOT)  “mEN 

CALL  WSCT(WEIGHT,NPTP1 ,GXX,GYY) 

CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 
CALL  PEAK(WORKX, IPNT, LPNT, PKVAL, PKLOC) 

IF  (SCTPK  .LT.  PKVAL)  THEN 
SCTPK  =  PKVAL 
IND-^-  PKLOC  -  IPNT  +  1 
SCTTAU(  ITRIAL)  =  FLOATdND  -  WlDTH/2) 
SCTTDT( ITRIAL)  =  TAUDOT 
END  IF 
END  IF 

IF  (AML)  THEN 

CALL  WAML(WEIGHT,NPTP1 ,GXYR,GXYI ,GXX,GYY) 
CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 
CALL  PEAK(V/ORKX,  IPNT, LPNT, PKVAL, PKLOC) 

IF  (AMLPK  .LT.  PKVAL)  THEN 
AMLPK  =  PKVAL 
IND  =  PKLOC  -  IPNT  +  1 
AMLTAU(  ITRIAL)  =  FLOATdND  -  WIDTH/2) 
AMLTDTd  TRIAL)  =  TAUDOT 
END  IF 
END  IF 
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IF  (MSC)  THEN 

CALL  WMSC(WEIGHT,NPTP1 ,GXYR,GXYI ,GXX,GYY) 

CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 

CALL  PEAK(WORKX, IPNT,LPNT,PKVAL,PKLOC) 

IF  (MSCPK  .LT.  PKVAL)  THEN 
MSCPK  =  PKVAL 
IND  =  PKLOC  -  IPNT  +  1 
MSCTAU(  ITRIAL)  =  FLOATdND  -  WIDTH/2) 

MSCTDT( ITRIAL)  =  TAUDOT 
END  IF 
END  IF 

IF  (OPT)  THEN 

CALL  WOPT(WEIGHT,NPTP1 ,FRQL,FRQU,FS,FTYPE,NSECT) 
CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 

CALL  PEAK(WORKX, I PNT,LPNT, PKVAL, PKLOC) 

IF  (OPTPK  .LT.  PKVAL)  THEN 
OPTPK  =  PKVAL 
IND  =  PKLOC  -  IPNT  +  1 
OPTTAU(  ITRIAL)  =  FLOATdND  -  WIDTH/2) 
OPTTDTd  TRIAL)  =  TAUDOT 
END  IF 
END  IF 

IF  (PH AT)  THEN 

CALL  WPHT(WEIGHT,NPTP1,GXYR,GXYI ) 

CALL  GCC(WORKX,WORKY,NFFT,WEIGHT,GXYR,GXYI ) 

CALL  PEAK(WORKX, IPNT, LPNT, PKVAL, PKLOC) 

IF  (PHTPK  .LT.  PKVAL)  THEN 
PHTPK  =  PKVAL 
IND  =  PKLOC  -  IPNT  +  1 
PHTTAU(  ITRIAL)  =  FLOATdND  -  WIDTH/2) 

PHTTDT( ITRIAL)  =  TAUDOT 
END  IF 
END  IF 

IF  (OUTPUT)  CLOSEd) 


C  '  ■  . 

C  GO  BACK  FOR  NEXT  DOPPLER  COMPENSATION 

C 

END  DO 
C 

C ************* ************************************** 

c*************************************************** 

c 

C  GO  BACK  FOR  NEXT  TRIAL 

C 

C 

END  DO 
C 

C************************************************************* 

c************************************************************* 

c 

C  OUTPUT  ESTIMATES  OF  TAU  &  TAUDOT 

C 
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IF  (SCO  THEN 

TNAM  =  'SCCTAU.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM, ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , SCCTAU , NTRL ) 

TNAM  =  'SCCTDT.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT • , TNAM, • REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , SCCTDT, NTRL ) 

END  IF 

IF  (SCOT)  THEN 

TNAM  =  'SCTTAU.'  //  SNRDB 

CALL  SGOPENd , 'WR I TE NOPROMPT ', TNAM, • REAL ', NTRL) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , SCTTAU , NTRL ) 

TNAM  =  'SCTTDT. '  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TN AM , ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , SCTTDT, NTRL ) 

END  IF 

IF  (AML)  THEN 

TNAM  =  'AMLTAU.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TN AM, ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , AMLTAU , NTRL ) 

TNAM  =  'AMLTDT. '  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM, ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , AMLTDT, NTRL ) 

END  IF 

IF  (MSC)  THEN 

TNAM  =  'MSCTAU.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ', TNAM ,' REAL ', NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , MSCTAU, NTRL ) 

TNAM  =  'MSCTDT. '  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM , ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , MSCTDT, NTRL ) 

END  I  F 

IF  (OPT)  THEN 

TNAM  =  'OPTTAU.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ', TNAM ,' RE AL ', NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , OPTTAU , NTRL ) 

TNAM  =  'OPTTDT. '  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM , ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , OPTTDT, NTRL ) 

END  IF 

IF  (PHAT)  THEN 

TNAM  =  'PHTTAU.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM , ' REAL ' , NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL  * , PHTTAU , NTRL ) 

TNAM  =  'PHTTDT.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ', TNAM ,' REAL ', NTRL ) 
CALL  SGTRAN ( 1 , ' WR I TE ' , ' REAL ' , PHTTDT, NTRL ) 

END  IF 

TNAM  =  'TDELAY.'  //  SNRDB 

CALL  SGOPEN ( 1 , ' WR I TE ' , ' NOPROMPT ' , TNAM, ' REAL ' , NTRL ) 
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CALL  SGTRAN ( 1 , *  WR I TE  * , *  REAL ' , TDEL AY , NTRL ) 

TNAM  =  'TDRATE.*  //  SNRDB 

CALL  SGOPEN ( 1 , *  WR 1 TE  * , *  NOPROMPT  * , TNAM , *  REAL ' , NTRL ) 

CALL  SGTRAN ( 1 , *  WR I TE ’ , ' REAL  * , TDRATE , NTRL ) 

CLOSE(I) 

C 

C 

C  GO  BACK  FOR  NEXT  SNR 
C 

END  DO 
C 

Q***it****************************************************************ititit 

C 

C  SAVE  SEED  VALUES  ON  DISK  IF  DESIRED 
C 

IF  (SVSEED)  THEN 

CALL  SSEED(SEED(1),SEED(7)) 

CALL  SGOPEN ( 98 , ' WR I TE  * , *  NOPROMPT ' , ' SEED . GCC ' , ' I NTEGER ' , 38 ) 
CALL  SGTRAN ( 98 , ' WR I TE  * , M  NTEGER ' , SEED , 38 ) 

END  IF 
C 

STOP 

END 
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Q************************************************** **********  *********** 


TIME  SHIFTED  SEQUENCE  GENERATOR 
VAX-11  FORTRAN  SOURCE  FILENAME: 
DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
REVISION  DATE 

00.0  JULY,  1983 


SHI  FT. FOR 

KANSAS  STATE  UNIVERSITY 
PROGRAMMER(S) 

K  SCARBROUGH 


********************************************************************** 


CALLING  SEQUENCE 

CALL  SH I FT( X 1 , X2, NPTS , NWT, DDOT, DELAY ) 

PURPOSE 

This  routine  generates  a  time  shifted  version  of  a 
sequence.  The  time  shift  is  linear  with  rate  DDOT 
for  the  length  of  the  sequence.  The  delay  rate 
can  be  changed  on  subsequent  calls  to  the  routine, 
however,  the  time  delay  is  computed  by  this  routine 
and  Is  continuous  between  calls. 

EMBEDDED  ROUTINE (S)  CALLED  BY  THIS  ROUTINE 

COEFF 

ARGUMENT(S)  REQUIRED  FROM  THE  CALLING  ROUTINE 

XI  -  Input  sequence  to  be  time-shifted 

NPTS  -  Number  of  points  in  Input  sequence 

NWT  -  Number  of  coefficients  for  'Shift-Filter' 

(NWT  is  assumed  to  be  odd,  and  .LE.  61) 

DDOT  -  Delay  rate 

DELAY  -  Initial  time  delay 

(used  only  on  first  call  to  subroutine) 

ARGUMENT! S)  SUPPLIED  TO  THE  CALLING  ROUTINE 

X2  -  Time-shifted  version  of  XI 

DELAY  -  Value  of  time  delay  after  NPTS 

*****«»»»*»»»»»»»*»»»***»*»*****»»»»»»******************* ************ 
C 

SUBROUTINE  SHIFT(X1 ,X2, NPTS, NWT, DDOT, DELAY) 

C 

REAL  X1(*),  X2(*),  XTMP(61),  WT(61),  PAST(61 ) 

LOGICAL  FIRST  /.TRUE./ 


¥ 
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SAVE  FIRST,  PAST,  TDEL 
C 

C************************************************************* 

c 

C  INITIALIZE  PARAMETERS 
C 

NWT2  =  ( NWT  -  1 )  /  2 
IF  (FIRST)  THEN 
FIRST  =  .FALSE. 

TDEL  =  DELAY 

CALL  RVINIT(PAST,NWT2,0.0) 

END  IF 
C 

DO  IPT  =  1,  NWT2 

JPT  =  IPT  +  NWT2 
XTMP(IPT)  =  PAST(IPT) 

XTMP(JPT)  =  XI (IPT) 

END  DO 
C 

C************************************************************* 

c 

C  GENERATE  TIME  SHIFTED  SEQUENCE 
C 

DO  IPT  =  1,  NPTS 
C 

0*************************************************** 

C 

C  CALCULATE  'TIME-SHIFT'  FILTER  COEFFICIENTS 

C 

CALL  COEFF(WT, NWT, TDEL) 

TDEL  =  TDEL  +  DDOT 
C 

C*************************************************** 

C 

C  IMPLEMENT  'TIME-SHIFT'  FILTER 

C 

KPT  =  IPT  +  NWT2 
XTMP(NWT)  =  XI (KPT) 

X2(IPT)  =  0.0 
DO  IWT  =  1 ,  NWT 

X2(IPT)  =  X2(IPT)  +  WT(IWT)  *  XTMP(IWT) 

END  DO 
C 

c 

C  SHIFT  VALUES  IN  XTMP  ARRAY 

C 

DO  I WT  =  1 ,  NWT  -  1 

XTMP( IWT)  =  XTMP( IWT  +  1 ) 

END  DO 
C 

END  DO 
C 

C 

C  SAVE  LAST  'NWT2'  VALUES  OF  XI  ARRAY 
C 
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DO  IPT  =  1,  NWT2 

JPT  =  IPT  +  NPTS  -  NWT2 
PAST(IPT)  =  XI (JPT)  . 

END  DO 
C 

Q *»«*#»»#»««»*»»»»**«****»««»*****************************»** * 


RETURN  TIME-SHIFTED  SEQUENCE  AND 
CURRENT  DELAY  VALUE  TO  MAIN  ROUTINE 

DELAY  =  TDEL 

RETURN 

END 
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SUBROUTINE  COEFF(V/T,NWT,TDEL) 

This  routine  computes  the  coefficients  for  a  filter 
which  generates  a  time  shifted  version  of  a  given 
sampled  sequence. 

NOTE  :  It  is  assumed  that  NWT  wi I  I  be  odd. 

«««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««« 

C 

SUBROUTINE  COEFF(WT,NWT,TDEL) 

C 

REAL  WT(*) 

C 

PI  =  3.141593 
NWT2  =  ( NWT  -  1 )  /  2 
C 

DO  IWT  =  1  ,  NWT 

JWT  =  NWT  -  IWT 

ARC  .=  (FL0AT(JWT  -  NWT2)  -  TDEL)  *  PI 
IF  (ARC. NE. 0.0)  THEN 

WT(IWT)  =  S INC  ARC)  /  ARC 

ELSE 

WT( IWT)  =1.0 
END  IF 
END  DO 
C 

RETURN 

END 
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